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A. Introduction 

Many physical problems involve calculating sums over paths. Each path could rep- 
resent one possible physical realization of an object such as a polymer, in which case the 
weight of the path is the probability of that configuration. The weights themselves could 
be imaginary as in the case of Feynman paths describing the amplitude for the propagation 
of a particle. Path integral calculations are now a standard tool of the theoretical physicist, 
with many excellent books devoted to the subject [|],0]. 

What happens to sums over paths in the presence of quenched disorder in the medium? 
Individual paths are no longer weighted simply by their length, but are influenced by the 
impurities along their route. The sum may be dominated by "optimal" paths pinned to the 
impurities; the optimal paths usually forming complex hierarchical structures. Physical 
examples are provided by the interface of the random bond Ising model in two dimensions, 
and by magnetic flux lines in superconductors. The actual value of the sum naturally 
depends on the particular realization of randomness and varies from sample to sample. I 
shall initially describe the problem in the context of the high temperature expansion for 
the random bond Ising model. Introducing the sums over paths for such a lattice model 
avoids the difficulties associated with short distance cutoffs. Furthermore, the Ising model 
is sufficiently well understood to make the nature of various approximations more evident. 

The high temperature correlation functions of the Ising model are dominated by the 
shortest paths connecting the spins. Such configurations, that exclude loops and overhangs, 
are referred to as directed paths. They dominate the asymptotic behavior of the sum over 
distances that are much longer than the correlation length. Most of the lectures are 
devoted to describing the statistical properties of sums over such directed paths. As in 
all multiplicative noise processes, the probability distribution for the sum is broad. Hence 
Monte Carlo simulations may not be an appropriate tool for numerical studies; failing to 
find typical members of the ensemble. Instead, we shall present a transfer matrix method 
that allows a numerical evaluation of the sum in polynomial time in the length of the path. 
The results indeed show that the sum has a broad probability distribution that resembles 
(but is not quite) log-normal. 

To obtain analytical information about this probability distribution we shall introduce 
the replica method for examining the moments. A brief review explains the relationship 
between the moments, the characteristic function, and cumulants. It can be shown easily 
that the one dimensional sum has a log-normal distribution. The moments of the sum 



over directed paths in two dimensions can be obtained by using a simple Bethe Ansatz. 
The implications and limitations of this approach are discussed. There is little analytical 
information in three and higher dimensions, but a variety of numerical results are available, 
mostly by taking advantage of a mapping to growing surfaces. 

The spin glass problem describes a mixture of ferromagnetic and antiferromagnetic 
bonds. The resulting sums in the high temperature expansion involve products over a 
random mixture of positive and negative factors. The calculation of moments is somewhat 
different from the case of purely positive random bonds. However, we shall demonstrate 
that the scaling behavior of the distribution is unchanged. A similar sum involving prod- 
ucts of random signs is encountered in calculating the probability of an electron tunneling 
under a random potential. In the strongly localized limit, it is again sufficient to focus on 
the interference of the forward-scattering (directed) paths. A magnetic field introduces 
random phases in the sum; while to describe the tunneling of an electron in the presence 
of spin-orbit scattering requires examining the evolution of a two component spinor and 
keeping track of products of random matrices. We shall argue that all these cases are in 
fact described by the same universal probability distribution which, however, does retain 
some remnant of the underlying symmetries of the original electronic Hamiltonian. 

Yet another class of directed paths has been introduced in the context of light scatter- 
ing in turbulent media. Assuming that inelastic scattering can be neglected, the intensity 
of the beam is left unchanged, and the evolution is unitary. Due to the constraint of uni- 
tarity the resulting directed paths are described by a probability distribution belonging to 
a new universality class. We shall introduce a discrete matrix model that explicitly takes 
care of the unitarity constraint. In this model, several properties of the resulting sum over 
paths can be calculated exactly. 

In the last sections of the course we shall go beyond the limitations of directed paths. 
In a uniform system the sum over all paths is calculated approximately, reproducing mean- 
field critical behavior. In two dimensions, the sum can be performed exactly for the Ising 
model. This leads to exact solutions for the pure Ising model, or other uniformly frus- 
trated two dimensional lattices. We shall use this method of exact summation to develop 
an integer algorithm for obtaining exact partition functions for two dimensional random 
lattices in polynomial time. Some results are described for the two dimensional spin glass, 
and contrasted with those obtained from Monte Carlo or transfer matrix methods. 

The results described in these notes were the outcome of many collaborations. In 
particular, I would like to express my thanks to E. Medina, L. Saul, Y. Shapir, and Y.-C. 



Zhang. I am grateful to F. David, P. Ginsparg, and J. Zinn-Justin for organizing the Les 
Houches summer school, and to the many students who helped me with correcting the 
lecture notes. The work at MIT was supported by the NSF through grants DMR-93-03667 
and PYI/DMR-89-58061. 

B. High Temperature Expansions for the Ising Model 

Consider a d dimensional hypercubic lattice of N sites. On each site there is an Ising 
variable o~i = ±1, and the spins interact through a Hamiltonian 

W = - 22 J iJ (T i (T J ■ W 

The symbol < ij > implies that the sum is restricted to the dN nearest neighbor bonds on 
the lattice. The bonds {Jij} are quenched random variables, independently chosen from a 
probability distribution p(J). For each realization of random bonds, the partition function 
is computed as 

z[j tj ] = J2 ex P H#0 = E II eKijaiaj > ( 2 ) 

{><} {ai}<ij> 

where the sums are over the 2^ possible configurations of spins, (3 = l/iksT) and K{j = 
(3 Jij. To obtain a high temperature expansion, it is more convenient to organize the 
partition function in powers of tanhi^j. Since (o~i<jj) 2 = 1, the Boltzmann factor for each 
bond can be written as 

e Kij _|_ e -Kij e Kij _ e -K i:j 



H ~ CTiCTj = COsh K i:j (1 + TijOiOj) , (3) 

z z 

where r^ = tanhKij is a good high temperature expansion parameter. Applying this 
transformation to each bond of the lattice results in 

W^^^^^EII^W;) . (4) 

Wi} Wi}{ij) 

where 

C dN = IjlcoshKij 

The term C dN is non-singular, and will be mostly ignored henceforth. The final 
product in eq. (H) generates 2 dN terms which can be represented diagrammatically by 



drawing a line connecting sites i and j for each factor of TijOiOj. Each site now obtains a 
factor of o\\ where < pi < 2d is the number of bonds emanating from i. Summing over 
the two possible values oi = ±1, gives a factor of 2 if pi is even and is pi is odd. Thus the 
only graphs that survive the sum have an even number of lines passing through each site. 
The resulting graphs are collections of closed paths Q on the lattice. The contribution of 
each graph is the product of r^ for the bonds making up the graph, resulting in 

Z[J lJ ] = 2 N xC dN J2\ II *■«] • ( 5 ) 

For a (i-dimensional hypercubic lattice the smallest closed graph is a square of 4 bonds 
and the next graph has 6 bonds. Thus, 



N w V^dN 



Z[Jij] = 2 JV x C 



1 + Yj TP1TP2TP3TP4 + 0(t 6 ) + ■ 



p 



(6) 



where the sum runs over the Nd(d — l)/2 plaquettes on the lattice and Tp a indicate the 
four bonds along each plaquette. A quench averaged free energy is now obtained as 



hi-Z , ,i ^-ft did — 1) A 

—— = ln2 + rflncoshK + ^- — ; -T 4 + --- , (7) 

where the over-lines indicate averages over the probability distribution p(J). 

The same method can be used to obtain expansions for various spin operators. For 
example the two spin correlation function is given by 

E, . .. KijUiOj ~7y HN 

(&m&n) = 2_^ y Cf m a n = —— 2_^ Cf m a n _[_[(!+ Tij^iVj) ■ (8) 

Wi} {at} (ij) 

The terms in the numerator involve an additional factor of a m a n . To get a finite value after 
summing over a m = ±1 and a n = ±1 we have to examine graphs with an odd number of 
bonds emanating from these external sites. After cancelling the common factors between 
the numerator and denominator, we obtain 



\<y m O n ) = ) r-^- . (9) 



Earn 



Ti 



Whereas the graphs in Q have an even number of bonds going through each site, those 
of Qmn have an odd number of bonds going through the external points m and n. This 
procedure can be generalized to multiple spin correlation functions. 



C. Characteristic Functions and Cumulants 

As equations (§) and (§) indicate, the partition function and correlation functions of 
the random system are themselves random quantities, dependent on all the bonds Kij. 
It may not be sufficient to just characterize the mean value of Z (or InZ), since the 
full information about these fluctuating quantities is only contained in their respective 
probability distributions p(Z) and p ((cr m a n )). It is thus important to learn to characterize 
and manipulate probability distributions, necessitating the short digression taken in this 
section to define various notations I shall use in describing such random quantities. 

Consider a continuous random variable x, whose outcome is a real number E. 

• The cumulative probability function (CPF) P(x), is the probability of an outcome with 
any value less than x, i.e. P(x) = prob.(i? < x). P(x) must be a monotonically increasing 
function of x, with P(— oo) = and P(+oo) = 1. 

• The probability density function (PDF), is defined by p(x) = dP(x)/dx. Hence, 
p(x)dx = prob.(x < E < x + dx). As a probability density, it is positive, and normalized 
such that 

/■OO 

P(oo) = / dxp(x) = 1 . (10) 



The expectation value of any function F(x), of the random variable is 

/"OO 

(F(x)) = / dxp(x)F(x) . (11) 



The function F(x) is itself a random variable, with an associated PDF of pp(f)df = 
prob.(/ < F(x) < f + df). There may be multiple solutions X{, to the equation F(x) = f, 
and 



PF(f)df = ^p(Xi)dXi => PF{f) = ^2p(xi) 



dx 



dF 



(12) 



The factors of \dx/dF\ are the Jacobians associated with the change of variables from x to 

F. 

• Moments of the PDF are expectation values for powers of the random variable. The 

n th moment is 

(x n ) = [ dxp{x)x n . (13) 



• The characteristic function, is the generator of moments of the distribution. It is simply 
the Fourier transform of the PDF, defined by 

p(k) = (e~ ikx ) = [ dxp(x) e~ ikx . (14) 



Moments of the distribution can be obtained by expanding p(k) in powers of k, 

«*> = (I) ^*") = E ^^V) ■ («) 

\n=0 ' / n=0 

• The cumulant generating function is the logarithm of the characteristic function. Its 
expansion generates the cumulants of the distribution defined through 

lnp(k) = jrt3l {x n )c . (16) 

n. 

Relations between moments and cumulants can be obtained by expanding the logarithm 
of p(k) in eq. (|H>1), and using 



ln ( 1 + e ) = E(- 1 ) n+1 - ■ ( 17 ) 



n=l 



The first four cumulants are called the mean, variance, skewness, and curtosis of the 
distribution respectively, and are obtained from the moments as 

(x) c = (x) , 

<z 2 ) c =<* 2 )-<*> 2 , 

(18) 
(x 3 ) c = (x 3 ) - 3 (x 2 ) (x) + 2 (xy , 

<x 4 ) c = (x 4 ) - 4 (x 3 ) (x) - 3 (x 2 ) 2 + 12 (x 2 ) (x) 2 - 6 (x) 4 . 

The cumulants are usually the most compact way of describing a PDF. An important 
theorem allows easy computation of moments in terms of the cumulants: Represent the n th 
cumulant graphically as a connected cluster of n points. The m th moment is then obtained 
by summing all clusters (connected or disconnected) of n points; the contribution of each 
cluster being the product of the connected cumulants that it represents. Using this result 
the first four moments are easily computed as 

(x) = (x) c , 






>=<* 2 > c +<*>; 



(19) 

<x 3 )=<x 3 ) c + 3<x 2 ) c (x) c + (x) 3 



' c 



2 



> = <x 4 ) c + 4 <x 3 ) c (x) c + 3 <x 2 » 6 <x 2 ) c (x) 2 + (x)\ 

7 



The normal (Gaussian) distribution describes a continuous real random variable x, 



with 



p(x) 



V2 



: exp 



7TCT Z 



(x ~ A)' 

2(7 2 



The corresponding characteristic function also has a Gaussian form, 



p(k) = 



dx 



1 



V2 



exp 



TIG" 



(x - xy 

2a 2 



ikx 



exp 



-ikX 



z 2 a 2 



(20) 



(21) 



Cumulants of the distribution can be identified from lnp(fe) = — ikX — k 2 a 2 /2, using 
eq. (|16|), as 



<*>c 



A 



<- 2 > ( 



a 



<- 3 > c = <- 4 > ( 







(22) 



The normal distribution is thus completely specified by its two first cumulants. This makes 
the computation of moments using the cluster expansion (eqs. (|19D) particularly simple, 

and 

(x)=X , 

X ' (23) 

(x 3 ) =3a 2 A + A 3 , 

< x 4 )=3a 4 + 6a 2 A 2 + A 4 , ••• . 

• The central limit theorem describes the probability distribution for a sum S = ^2 i=1 Xi 
over a large number of random variables. Cumulants of the sum are given by, 

N N 

<S>c = £>>c . <^> c = E(^)c . ■■■ ■ (24) 



i=i 



'■■3 



If the random variables are independent, p(x) = [jp^i;), and ps(k) = YiPi(k). The cross- 
cumulants in eq. (|2^ ) vanish, and the n th cumulant of S is simply the sum of the individual 
cumulants, (<S n ) c = J2i=i ( x7 i)c- When all the iV random variables are independently 
taken from the same distribution p(x), this implies (<S n ) c = N (x n ) c . For large values 
of N, the average value of the sum is proportional to A^ while fluctuations around the 
mean, as measured by the standard deviation, grow only as viV. The random variable 
y = (S — N (x) c )/vN, has zero mean, and cumulants that scale as (y m ) c oc N l ~ m / 2 . As 
A^ — > oo, only the second cumulant survives and the PDF for y converges to the normal 
distribution, 



lim p J y 

N—i-oo 






N(x) c 



1 



VM?l 



: exp 



y 



2(x 2 ) ( 



(25) 



The convergence of the PDF for the sum of many random variables to a normal distribution 
is a most important result in the context of statistical mechanics where such sums are 
frequently encountered. The central limit theorem proves a more general form of this 
result: It is not necessary for the random variables to be independent, as the condition 



>A f 



Y^i ... i (xi t ■ ■ ■ Xi m ) c "C 0(N m > ), is sufficient for the validity of eq. 



D. The One Dimensional Chain 

The graphical method provides a rapid way of solving the Ising model at zero field 
in d = 1. We can compare and contrast the solutions on chains with open and closed 
(periodic) boundary conditions. 

1. An open chain of iV sites has N — 1 bonds. It is impossible to draw any closed graphs 
on such a lattice, and hence 



N-l 



Z = 2 N f[ coshK Q x 1 



a = l 



InZ 

~n~ 



In [2 cosh K] 



ln[coshiir] 

N 



(26) 



where K a = K aa+ i. There is also only one graph that contributes to the two point 
correlation function, 



\o- m o- n } = 2_^ y o- m a r 



n-l 

n 



TV, 



(27) 



2. A closed chain has the same number of sites and bonds, N. It is now possible to draw 
a closed graph that circles the whole chain, and 



A f 



Z = 2 N Y[ coshK c 



va=l 



IV 



i+ n 



TV, 



vo=l 



(28) 



InZ 



In [2 cosh K] + — r 



N 



The difference between the quenched free energies of closed and open chains is a surface 
term of the order of 1/N, and an exponential decay reflecting the interaction between 
edges, both vanishing in the thermodynamic limit oi N —>■ oo. The correlation function 
can again be calculated from eq. (H). There are two paths connecting the points m and n, 
along the two possible directions on the chain, giving 



(0'mO'n) 



n 



n-l 
a=?7t a 



+ n 



771 — 1 



N 



1+ IEUra 



(29) 



Since the partition function of the open chain is the sum of N— 1 independent variables 

InZ 



N 



^N-i In cosh K a 
hi 2 + > ^7 



(30) 



we can use the central limit theorem to conclude that as N — > oo the probability distribu- 
tion p(lnZ/N), is a gaussian with mean 



ln Z = N (ln 2 + ln cosh K) 



and variance 



(ln Zf c = (ln Z) 2 - (ln Z) = AT (In cosh X) 2 . 



(31) 



(32) 



(Note that I have ignored the small difference between N and N — 1 in the thermodynamic 
limit.) Similarly, for the correlation function of two points separated by a distance t, we 

have 

Et-i 
lnr a . (33) 

a=0 

As long as the random variables on the bonds are independently distributed, the cumulants 
of ln (<jQ<7t) are given by, 



ln(o"ocr t ) = t IntanhAT 



(ln(a-ocx t )) 2 =t (IntanhK)l 



(34) 



k {]n(er (Tt))c = t (lntanhK)! 



In the following sections we shall try to obtain similar information about probability 
distribution functions for the partition and correlation functions in higher dimensions. To 
do this we shall employ the replica method for calculating the moments of the distribution. 
For example, the cumulants of the free energy are given by 



exp 



n 



IV 



nlnZ+— (In Z)2 + ... + — (In Z)§ + 



pi 



(35) 



where we have taken advantage of eq. ([H|), replacing (—ik) with n. Usually, the moments 
on the left hand side of the above equation are known only for integer n, while the evalu- 
ation of the cumulants on the right hand side relies on an expansion around n = 0. This 
is one of the difficulties associated with the problem of deducing a probability distribution 
p(x), from the knowledge of its moments x n . There is in fact a rigorous theorem that 



10 



the probability distribution cannot be uniquely inferred if its n th moment increases faster 
than n! 0. Most of the distributions of interest to us (such as the above log-normal) do 
not satisfy this condition! Similar problems are encountered in the replica studies of spin 
glasses 0] . It turns out that many of the difficulties associated with a rigorous inversion 
are related to the tail of the distribution. Most of the information of interest to us is con- 
tained in the "bulk" of the distribution which is easier to investigate. Rather than taking 
a rigorous approach to the problem, we shall illustrate the difficulties and their resolution 
by examining the one dimensional case in detail since it actually presents the worst case 
scenario for the inverting of moments. 

We used the central limit theorem to deduce that the probability distribution for 
((Jo&t) is log-normal. Its moments are computed from, 



t-i 



exp 



a=0 



p 



(36) 



Let us consider a binary distribution in which r takes two positive values of T\ and r 2 > T\ 

with equal probability. Then 

^=lL±ll_ f (37) 



and the generating function for the cumulants of the correlation function is 

ln^=nlnr 1 + ln' 1 + (r2/ri)n ' 



n 



n In n + In 



=n In n + In 



l + l + nln(r 2 /ri) + n 2 /21n 2 (t 2 /ti) + 



^H? + 7"; l + 



(38) 



-■»- + HS + t* , £ 



+ • 



=nln( v /rir 2 ) + — In 

■ n 



+ • 



Combining eqs. ( |36|) and (p8|), the cumulants of the correlation function are given by 



ln(a a t ) =*ln(0-ir 2 ) 



]n(a cr t )c =7 ln (^/ r l) 



(39) 
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While it is true that ln(a"ocr t ) is normally distributed for large t, we should be careful 
about the order of limits in terminating the power series in the exponent at the second 
order. If we do so, from 



(aoa t ) n ~ exp [t (nai + n 2 a 2 )] , (40) 

we should not infer anything about the high moments (n — > oo) and the tail of the dis- 
tribution. Otherwise (since a 2 > 0), we would conclude that sufficiently large moments of 
(o"ocr t ) diverge with separation; a clearly false conclusion as (o"oa t ) is bounded by unity! 
The exact result is that 

nt 

lim (*o<rt) n = if > ( 41 ) 

n — >oo Z 

i.e., the high moments are almost entirely dominated by the one exceptional sample in 
which all bonds are equal to t 2 . We can summarize the situation as follows: The "bulk" of 
the probability distribution for ln(o"ocr t ) is described by the small moments (n — ► 0), while 
the tail of the distribution is governed by the large moments (n — ► oo). We should have 
a clear idea of the crossover point n* in applying the replica method. For the above one 
dimensional example, an estimate of n* is given by the ratio of the successive terms in the 
expansion, i.e. 

"* = w 1 / x • (42) 

ln(r 2 /ri) 

Note that as rijr\ becomes large, n* decreases, possibly becoming smaller than unity. 
This does not imply that we should conclude that ln(<7o0t) is not normally distributed, 
just that the tail of the distribution is more prominent. Failure to appreciate this point is 
the source of some misunderstandings on the use of the replica method |J . 

Clearly, it is possible to come up with many different microscopic distributions p(r), 
which result in the same first two cumulants in eqs. fl39[), but different higher cumulants. 
All these cases lead to the same universal bulk probability distribution for ln(a" cr t ) at 
large t, but very different tails. Thus the non-uniqueness of the overall probability in this 
example has to do with the rather uninteresting (and nonuniversal) behavior of the tail 
of the distribution. The correct interpretation of eqs. (|39| ) is that the mean value for the 
logarithm of the correlation function grows linearly with the separation t. In analogy with 
pure systems, we can regard the coefficient of this decay as the inverse correlation length, 
i.e. £ _1 = —ha.y/TiT2- However, due to randomness in the medium, correlations have 
different decays between different realizations (and between different points in the same 
realization). The variations in this "inverse correlation length" are scale dependent and 
fall off as 1/y/i. In the next sections we shall attempt to generalize these results to higher 
dimensions. 

12 



E. Directed Paths and the Transfer Matrix 

Calculation of the correlation function in higher dimensions is complicated by the 
presence of an exponentially large number of paths connecting any pair of points. On 
physical grounds we expect the high temperature phase to be disordered, with correlations 
that decay exponentially as a function of the separation t. The essence of this exponential 
decay is captured by the lowest order terms in the high temperature expansion. The first 
term in the series comes from the shortest path connecting the two points. Actually, along a 
generic direction on a hypercubic lattice there are many paths that have the same shortest 
length. (In two dimensions, the length of the shortest path connecting (0, 0) to (t, x) is the 
"Manhattan" distance \t\ + \x\.) The number of paths grows from a minimum of 1 along a 
lattice direction to a maximum of d per step along the diagonal. (The number of paths on 
the square lattice is (t + x)l / (tlxl) .) Thus the decay of correlations depends on orientation, 
a consequence of the anisotropy of the hypercubic lattice. (Note that this anisotropy is 
absent at distances less than the correlation length. We thus don't have to worry about 
anisotropy effects in discretizing critical (massless) theories on a lattice.) 

In a uniform system these shortest paths are sufficient to capture the essence of cor- 
relation functions of the high temperature phase: An exponential decay with separation 
which is generic to all spin systems. As temperature is reduced, more complicated paths 
(e.g. with loops and overhangs) start contributing to the sum. Although the contribution 
of these paths decays exponentially in their length, their number grows exponentially. Ul- 
timately at the critical point this "entropic" increase in the number of paths overcomes the 
"energetic" decrease due to the factors of r < 1, and paths of all length become important 
below T c . However, throughout the high temperature phase, it is possible to examine the 
paths at the coarse grained scale where no loops and overhangs are present. The scale 
of such structures is roughly the correlation length £, and if we use £ as the unit of a 
coarse-grained lattice, the paths contributing to the correlation function are directed. 

Let me define "directed paths" more carefully: Between any pair of points on the 
lattice we can draw an imaginary line which I shall refer to as the "time" axis t. Transverse 
directions (perpendicular to the t axis) are indicated by x. Directed paths are similar to 
the worldlines of a particle x(t) in time; they exclude any path from the initial to the final 
point that has steps opposite to the main time direction. The question of the validity 
of this approximation, and the importance of the neglected loops continually comes up. 
This is possibly because it is more common to think about the vicinity of the critical 
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point where loops of all sizes are present and equally important. Away from the critical 
point, we have to distinguish between properties at scales smaller and larger than the 
correlation length £; there is no other length scale (except for the lattice spacing) in the 
problem. Limiting the sum for the correlation function to directed paths is only useful for 
separations i > ^. Loops, overhangs, and additional structures occur up to size £ (the 
only length scale present) and can be removed by coarse graining such that the lattice 
spacing is larger than or equal to £. This is automatically the case in a high temperature 
expansion since £ is initially less than a lattice spacing. By the same argument, we may 
also neglect the closed loops (vacuum bubbles) generated by the denominator of eq. (§). 

In this section I shall demonstrate how sums over directed paths in the uniform system 
can be calculated exactly by transfer matrix methods. The method can also be generalized 
to random systems, providing an algorithm for summing all paths in polynomial time. For 
ease of visualization, I shall demonstrate this method in two dimensions; the results are 
easily generalizable to higher dimensions. Also to emphasize the general features of the 
transfer matrix method, we shall compare and contrast the behavior of correlations along 
the axis and the diagonal of the square lattice. 

To calculate the correlation function (o"o,o<3"o,t), on a non-random square lattice, we 
shall focus on directed paths oriented along the main axis of the square. These paths 
are specified by a set of transverse coordinates (xq, Xi, X2, ■ ■ ■ , Xt), with xq = xt = 0. Of 
course, there is only one shortest path with all Xi equal to zero, but we would like to 
explore the corrections due to longer directed paths. Consider the set of quantities 

(x, t\W\0, 0) = sum over paths from (0, 0) to (x, t) = W(x, t) . (43) 

The calculation of W(x, t) is easily accomplished by taking advantage of its Markovian 
property: Each step of a path proceeds from its last location and is independent of the 
previous steps. Hence W can be calculated recursively from, 

W(x, t+l) = r \W(x, t)+r (W(x - 1, t) + W(x + 1, t)) + 0(r 2 )] 

v- (44) 

= Y, xl {AT\x')W{x',t) , 

where we have introduced a transfer matrix, 

(x\T\x')=t5 X!X/ +t 2 (S x ^ +1 +5 x ^- 1 ) + 0(t 3 ) . (45) 
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If we treat the values of W at a particular t as a vector, eq. (| 
W{t) = TW(t -l) = --- = T t W(0) 
starting from 



can be iterated as, 



(46) 



W(0) 



T 
1 



W 



(47) 



The calculations are simplified by diagonalizing the matrix T, using the Fourier basis 
(x\q) = e iq - x /\fN, as 

T(q) = r (1 + 2rcosg + •••)= rexp 



2TU-1T + 



(48) 



In this basis, W is calculated as 

w(x,t) = (x|r*|o) = Y^(x\Q) T (qY(q\o) 



T t e 2rt 



f dq r . 2 -, 

/ — - exp \iqx - q rt H 

J l-K 



(49) 



exp[t(lnr + 2r + C(r 2 ))] 



exp 



x 



Art 



The result is proportional to a gaussian form in x of width \/1rt. The exponential decay 
with ^ _1 = ln(l/r) — 2r + 0(r 2 ) at x = is accompanied by a subleading 1/y/i. 

The corresponding calculation of paths along the diagonal, contributing to (oo^Co,*) , 
is even simpler. (Note that the t and a; axes are rotated by 45° compared to the previous 
example.) At each step the path may proceed up or down, leading to the recursion relation 

W(x,t+l)=T(W(x-l,t) + W(x + l,t)) = J2( x \ T \ x ') W ( x '^) > ( 50 ) 

x' 

with the transfer matrix 

(x\T\x') = r (S X)X >+i + S X)X f-i) =4> T(q) = 2rcosq . (51) 

The calculation of W proceeds as before, 

W(x,t) = (x|T*|0) = V (x\q)T( q y (q\0) 

= /£(2-) i (cos,)V- j (52) 



(2t)* x 



1 
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exp 



x 

Yt 
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where the final result is obtained by a saddle point evaluation of the integral, essentially 
replacing cos* q with exp (—q 2 t/2) . 

The similarity between eqs. (|49|) and (l52| ) is apparent. Note that in both cases the 
leading exponential decay is determined by T(q = 0), i.e. 

W(0,t) « ALx = T(q = 0Y . (53) 

This is an example of the dominance of the largest eigenvalue in the product of a large 
number of matrices. There is a corresponding ground state dominance in the evolution of 
quantum systems. The similarities become further apparent by taking the continuum limit 
of the recursion relations, which are obtained by regarding W(x,t) as a smooth function, 
and expanding in the derivatives. From eq. (HI), we obtain 



W+^ + ...= tW + t*(2W + ^ + ...) , (54) 



while eq. flSCf ) leads to 



dW d 2 W 

W+ _ + ... = 2tW + t — + ... . ( 55 ) 



For large t, the function W decays slowly for adjacent points in the x direction, and it is 
justified to only consider the lowest order derivatives with respect to x. The decay factor 
along the t direction is, however, quite big and we shall keep track of all derivatives in this 
direction, leading to 

e dt W = t exp [2r + rdl + ■■■]}¥ , (56) 



and 



e dt W = 2rexp 



\ol + 



W , (57) 



respectively. Both equations can be rearranged (and generalized in higher dimensions) into 

the differential form, 

dW W 

where £(#) and v(9) are the orientation dependent correlation length and dispersion coef- 
ficient. Eq. ( J58D can be regarded as a diffusion equation in the presence of a sink, or an 
imaginary-time Schrodinger equation. 
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It is of course quite easy to integrate this linear equation to reproduce the results in 
eqs. (^9]) and (|52j). However, it is also possible]]]] to express the solution in the form of a 
continuous path integral. The solution is trivial in Fourier space, 

dW(q) 



dt 



-(C +vq)w 



W(q,t + At) = e-G 1+uq2)At W(q,t) 



(59) 



while in real space, 

W(x,t + At) 



dq 

2^ 



e^e-K 1+ucl ^ At W(q,t) 



—e iqx e 

2% 


'^~ 1 -^ M fdx t 


s -iQxt W ( 


%t,t) 


dxt exp 


At (x — xt) 
£ 4vAt 


W(x t ,t) 


dxt exp 


At At fxt+At -xA 
£ 4i/ V At J 


W( 



(60) 



which is just a continuum version of eqs. (H) and (|50|). We can subdivide the interval 
(0, t) into N subintervals of length At = t/N. In the limit of iV — » oo, recursion of eq. 
gives 



W(x,t) 



(x,t) 



Vx{t') exp 



(0.0) 



dt' 



1 



X 



£{0) Au{9) 



(61) 



where x = dx/dt', and the integration is over all functions x(t'). 

It is instructive to compare the above path integral with the partition function of a 
string stretched between (0,0) and (x,t), 

.(x,t) 



Z(x,t) 



Vx{t') exp 



13a I dtW^ + x 2 

(0,0) Jo 

(x,t) r rt 

Vx{t') exp 

(0,0) 



r*'0+ir 2+ -) 



(62) 



where a is the line tension. Whereas for the string £ -1 = (2z/) _1 = (3a, in general due to 
the anisotropy of the lattice these quantities need not be equal. By matching solutions at 
nearby angles of 9 and 9 + d9, it is possible to obtain a relation between £ -1 (#) and u~ 1 (9). 
(For a similar relation in the context of interfaces of Ising models, see ref. ||.) However, 
£(#) -1 calculated from the shortest paths only is singular along the axis 9 = 0. This is why 
to calculate the parameter v along this direction it is necessary to include longer directed 
paths. 
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F. Moments of the Correlation Function 

We now return to the correlation functions in the presence of random bonds. In the 
high temperature limit, we can still set 

t 

P i=l 

where the sum is over all the diagonally oriented directed paths P from (0, 0) to (x, t), and 
the Tpi denote factors of tanh K encountered for the random bonds along each path. The 
T% are random variables, independently chosen for each bond. We shall assume that the 
probability distribution p(r) is narrowly distributed around a mean value r with width 
a. Clearly, W(x, t) is itself a random variable and we would like to find its probability 
distribution. Rather than directly calculating p(W), we shall first examine its moments 



W n . 

Calculation of the first moment is trivial: Each factor of r% occurs at most once in 
eq. (^3|), and hence after averaging, 

W{x~t) = (ao,oCT x ,t) = J2 Tt ■ ( 64 ) 

p 

This is precisely the sum encountered in a non-random system, with r replacing r. For 
example, along the square diagonal, 



W(x, t) « (2r) f x -^= exp 



9 

ar 
2t 



(65) 



and in general, in the continuum limit, 



dW W d 2 W 



dt £ dx 2 

For the calculation of the second moment we need to evaluate 

t 



WW=J2H T ^^ ■ ( 67 ) 

p,p' i=i 

For a particular i, there are two possible averages depending on whether or not the two 
paths cross the same bond, 

J> 2 if Pi^P'i / rfi N 

T Pl T p/i = <— i- . (68) 

[ T A if Pi = Pi 
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Since t 2 > T 2 there is an additional weight for paths that intersect compared to those 
that don't. This can be regarded as an attraction between the two paths, represented by 
a Boltzmann weight, 



U= -* 



^2 y 2 T 2 



(69) 



Including the attraction, the recursion relation for WW is, 

W 2 (xi,x 2 ,t) = W(x u t)W(x 2 ,t) =J2 , l {xiX2\T 2 \x' 1 x' 2 )W{x' 1 ,x' 2 ,t-l) , (70) 

x \ x 2 

with the two body transfer matrix 



Xi 



;ix 2 |T 2 |x / 1 X2) = r 2 (5 Xl , x > i+1 + S Xl ,x' x -i) {$x 2 ,x> 2 +i + ^x 2 ,x' 2 -i) (l + (U - l)S XljX2 5 x ^ x ' 2 ) . 

(71) 
The significance of the attraction in eq. fl69|) is as follows: In the random system the 
paths prefer to pass through regions with particularly favorable values of r. After per- 
forming the quench averaging the paths go through a uniform medium. The tendency for 
the original paths to bunch up through favorable spots is instead mimicked by the uniform 
attraction which tends to bundle together the paths representing the higher moments. 
In the continuum limit, eq. ( |7D| ) goes over to a differential equation of the form, 

dW 2 (x x ,x 2 ,t) 2W 2 d 2 W 2 2 W 2 



at 



+ v—^ + v— 1 ^+u5{x 1 -x 2 )W 2 = -H 2 W 2 , (72) 



dx\ 



£ ' dx\ 

with u w a 2 /r 2 . Alternatively, we could have obtained eq. ([72] ) from the continuum 
version of the path integral, 



(Xl,t) 



W 2 (x ± ,x 2 ,t) = I Vx x (t') Vx 2 (t')exp 

'(0.0) J (0.0) 



(X2,t) 



exp 



dt' 






I dt'uS (xi(t') — x 2 (t')) 
Jo 


exp 


[s:<-\-m 



(73) 



Formally integrating eq. ([72] ) yields W 2 oc exp(— tH 2 ), which can be evaluated in the basis 
of eigenvalues of H 2 as 



W 2 (xi,x 2 ,t) = (xiX2|T 2 *|00) = ^(xix 2 |m)e- em *(m|00) 



-e t 



t — > OO 



(74) 



where {Em} are the eigenenergies of H 2 , regarded as a quantum Hamiltonian. The expo- 
nential growth of W 2 for t — ► oo is dominated by the ground state Eq. 
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The two body Hamiltonian depends only on the relative separation of the two particles. 
After transforming to the center of mass coordinates, 



r —x\ — X2 
R=(x 1 + x 2 )/2 
the Hamiltonian reads, 



dl + dl = \dl + <20, 



2 

r i 



H2 = j- U -d 2 R - 2ud 2 r - uS(r) 



(75) 



(76) 



The relative coordinate describes a particle in a delta-function potential, which has a 
ground state wavefunction 

Mr, R) oc e~ K \ r \ . (77) 

The value of k is obtained by integrating H 2 ^q from — e to e, and requiring the discontinuity 
in the logarithmic derivative of i/jq to match the strength of the potential; hence 



-2v{— k — k) = u 



u 



a 



K 



4v 2t : 



The ground state energy of this two particle system is 



The inequality, 



W 2 (t) =exp 



£q = +7 — 2l/K z 



it uH 



u 



« %!> 



T + 



uH 



W{t) l exp — - > W(t) 



(78) 



(79) 



(80) 



implies that the probability distribution for W(t) is quite broad, and becomes progressively 

wider distributed as t — > oo. 

Higher moments of the sum are obtained from 

t 

W^= E U T ^--- T ^ ■ ( 81 ) 

P 1 ,-,P n i=l 

At a particular "time" slice there may or may not be intersections amongst the paths. Let 

us assume that r is Gaussian distributed with a mean r, and a narrow width a; then, 



— f dx x m 

T m PS / ; exp 



V2 



na" 



(x -rf 
2a 2 



(set x = t + e and expand in e) 



de 



V2no^ 



_ m + mT m-i e + m(m 1) _ m _ 2(2 + 



exp 



2<r 



-m | m{m-l) _ 2 2 
t H r <7 + 



m(m — 1) a 2 



T 2 



(82) 



r m exp 



m(m — 1) 



-?i 
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Since there are m(m — l)/2 possibilities for pairing m particles, the above result represents 
the Boltzmann factor for a pairwise attraction of u for particles in contact. Since r is 
bounded by unity, the approximations leading to eq. (|82]) must break down for sufficiently 
large m. This implies the presence of three and higher body interactions. Such interactions 
are usually of higher order and can be safely ignored. For a discussion of these higher order 
interactions in a similar context see ref. 0. 

The continuum version of the resulting path integral is 



W n (x ir ■ ■ ,x„,t) =W(xi,t) ■ ■ -W(x n ,t) = / Vx x {t')---Vx n {t') 

J(0,0,—,0) 



exp 



and evolves according to 



nt 



t- dt ' EiS-fE^'')-^')) 



(83) 



dt 



nW n -A d 2 W n u sr^ , , , v r witt tj w 

+ v Z_^ ~f)~T~ + 9 2^ S ( x «( r ) ~ x p( T )) W n = -H n W n 



c 



1 dx l 

Q = l " 



a^/3 



(84) 



The asymptotic behavior of W n at large t is controlled by the ground state of H n . The 
corresponding wavefunction is obtained by a simple Bethe Ansatz Q, which generalizes 
eq. ([TFD to 



ipo(x 



w • • 



x n ) oc exp 



V — V 



(/ 



with k = — 
Av 



(85) 



For each ordering of particles on the line the wave function can be expanded as 
■i/>o oc exp [kqIq], with the "momenta" K a getting permuted for different orderings. For 
example, if x\ < £2 < ■ ■ ■ < x n , the momenta are 



K a = k [2a — (n + 1)] 



(86) 



forming a so called n-string. The kinetic energy is proportional to 

s=y n \2a-(n + i)} 2 = y n r(™ + i) 2 - 4 «(™ + i)+ 4 « 2 l 



=n(n + l) 2 -4(n + l) 



n(n + l) 4n(n+l)(2n+l) 



+ 



=n(n + 1) 



-(n+l) + 



2(2n + l) 



6 

n(n + l)(n — 1) 



(87) 
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leading to the ground state energy 

n o 

£ o = T-vJ2 K « = J-^- n ( n - 1 ) • ( 88 ) 

Thus the asymptotic behavior of moments of the sum has the form 



lim W n (t) = exp 



nt n(n — 1)uk t 



t 3 



n(n 2 — l)u 2 t 



W{t) n exp[ v ^J ) . (89) 



G. The Probability Distribution in Two Dimensions 

It is tempting to use eq. ( |89D in conjunction with 

(lnW)2 + --- + 

pi 



]im hi \W n (t)) =nlnW + — (In W) 2 C + ••• + — (In W) p c + ■ ■ ■ , (90) 



to read off the cumulants for the probability distribution for InVF. The key point is 

the absence of the n 2 term and the presence of the n 3 t factor in the exponent of eq. (|89|), 

suggesting a third cumulant, and hence fluctuations in In W that grow as t 1 ' 3 0. However, 

as discussed before, there are subtleties in trying to deduce a probability distribution from 

the knowledge of its moments which we need to consider first. Since W(t) is bounded 

by unity, eq. (J89Q cannot be valid for arbitrarily large n. Our first task is to identify the 

crossover point n* beyond which this result is no longer correct. 

Eq. fl88D is obtained for the ground state of n particles subject to a two body interaction 

in the continuum limit. A simple argument can be used to understand the origin of the n 3 

term in the energy, as well as the limitations of the continuum approach. Let us assume 

that the n particles form a bound state of size R. For large n, the energy of such a state 

can be estimated as 

n vn un 2 .„„. 

e ~l + w~ ■ (91) 

A variational estimate is obtained by minimizing the above expression with respect to R, 
resulting in R oc v/{un) and e oc u n /v. The size of the bound state decreases with 
increasing n, and the continuum approximation breaks down when it becomes of the order 
of the lattice spacing for n* oc v/u ~ T 2 /a 2 . For n ^> n* all the paths collapse together 
and 



lim W n (t) ~ (2r n ) . (92) 

n — >oo 
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This asymptotic behavior is non-universal and depends on the extreme values of the local 
probability distribution for r. Depending on the choice of parameters, n* can be large 
or small. However, as discussed in the context of the one dimensional problem, its value 
controls only the relative importance of the tail and the bulk of the probability distribution 
for In W. The behavior of the bulk of the distribution is expected to be universal. The 
crossover at n* is explicitly demonstrated in a simpler model in ref. ||10|| . 

Another important consideration is the order of limits. Eq. ( |B9"1) is obtained by taking 
the t — > oo limit at fixed n, while the cumulant series in eq. ([JDp relies on an expansion 
around n — > for fixed t. The two limits do not commute. In fact, we would naively deduce 
from eq. (|5pjj that only the third cumulant of In W is non-zero. This is incorrect as it is 
impossible to have a probability distribution with only a third cumulant [IT | . The correct 



procedure |L(J is to assume that the singular behavior associated with n — > and t — > oo 
is described by a scaling function of the form g s (nt u ). (This is similar to the singularity 
of the free energy at a critical point with n — 1 ' w playing the role of a correlation length.) 



If t — > oo at fixed n, extensivity of the free energy of the n particle system forces In W n {t) 
to be proportional to t. At the other limit of n —>■ at fixed t, the result is a power series 
in n, i.e. 

ant + pn x l^t for t — > oo at fixed n 



\nW n {t) = ant + g 8 (nt w ) = {""" ^ > JIL " „ ^ * ^ ^ tt " "*™ '" .(93) 

w ysv ; \ant + g x nt u + g 2 (nt") H f °r fixed t as n -> v ; 

Note that I have included a non-singular term, ant. Similar considerations have been put 



forward in ref. [12]. Comparison with eq. (|89|) gives oo = 1/3, and we can read off the 



cumulants of In W as 

]nW(t) =at + g x t 1/3 

\nW 2 (t) c =2g 2 t 2/3 

(94) 

{\nWP(t) c =p\g p tP/ 3 

The existence of t 1 ' 3 corrections to the quench averaged value of lnVF(t) was first 
suggested by Bouchaud and Orland ]13[] and has been numerically verified |lj§. The t 2 ! 3 
growth of the variance of the probability distribution was obtained by Huse and Henley 
||15|| in the context of interfaces of Ising models at zero temperature where an optimal path 
dominates the sum. The results remain valid at finite temperatures [ |16| . Simulations are 

23 



performed by implementing the transfer matrix method numerically. For example, along 
the diagonal of the square lattice, the recursion relation 

W(x,t+l)=T x ,t t -W(x-l,t)+T Xtti+ W(x + l,t) , (95) 

is iterated starting from W(x,0) = 5 Xt Q. The random numbers r Xtt ,a are generated as the 
iteration proceeds. The memory requirement (the arrays W(x)) depend on the final length 
t) each update requires t operations, and the total execution time grows as t 2 . Thus for a 
given realization of randomness, exact results are obtained in polynomial time. Of course 
the results have to be averaged over many realizations of randomness. The typical values of 
t used in the transfer matrix simulations range from 10 3 to 10 4 , with 10 2 to 10 3 realizations. 
Calculating higher cumulants becomes progressively more difficult. The existence of the 
third cumulant was verified by Halpin-Healy [17|. A fourth cumulant, growing as t 4 ' 3 
was observed by Kim et al. JIB] . Starting from the replica result, Zhang [jn| proposed 
an analytical form, p(lnW,t) ~ exp(— a| InW — In W\^' 2 /t 1 ' 2 ). While this form captures 
the correct scaling of free energy fluctuations, it is symmetric about the average value 
precluding the observed finite third cumulant. This deficiency was remedied by Crisanti 
et al. |2(J who generalized the above probability to one with different coefficients a± on 
the two sides of the mean value. 

So far, we focused on the asymptotic behavior of W(x,t) at large t, ignoring the 
dependence on the transverse coordinate. For the pure problem, the dependence of W 
on the transverse coordinate is a Gaussian, centered at the origin, with a width that 
grows as t 1 ' 2 . The full dependence is obtained in the pure problem by including the band 
of eigenvalues with energies close to the ground state. Unfortunately, determining the 
appropriate eigenvalues for the interacting problem is rather difficult. In addition to the 
eigenvalues obtained by simply multiplying eq. ( |5Bj ) by exp [iq [x\ + ■ ■ ■ + x n )}, there are 
other states with broken replica symmetry [21]]. A treatment by Bouchaud and Orland 
jnj includes some of the effects of such excitations but is not fully rigorous. It does 
predict that the extent of transverse fluctuations grows as ft with ( = 2/3 as observed 
numerically [|T5| , |T6]1 ■ There is in fact a relation between the exponents £ and uj which follows 
from simple physical considerations ||15| : By analogy with a string, the energy to stretch 
a path by a distance x grows as x 2 /t. The path wanders away from the origin, only if 
the cost of this stretching can be made up by favorable configurations of bonds. Since the 
typical fluctuations in (free) energy at scale t grow as t w , we have 

2 

— oc t u => u = 2C - 1 . (96) 
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This relation remains valid in higher dimensions and has been verified in many numerical 
simulations. The first (indirect) proof of a; = 1/3 was based on a replica analysis of a 
problem with many interacting paths ||22|| . It was soon followed by a more direct proof 



|23fl based on a completely different approach: the Cole-Hopf transformation described in 



next section. 



H. Higher Dimensions 



The approach described in the previous sections is easily generalized to higher dimen- 
sions. The directed path in d = D + 1 is described by x(t), where x is a D dimensional 
vector. The recursion relation of eq. (031) is generalized to 



W(x,t + l) = J2r £ -eaW(x-e l ,t) , (97) 

where e*j are unit vectors. The recursion relation is easily iterated on a computer, but 

the memory requirement and execution time now grow as t D and t D+1 respectively. The 
continuum limit of this recursion relation is 

a ^m = .^ + ^ w+ ^ t)w , (98 ) 

where fu.(x,t) represents the fluctuations of r(x,t) around its average. Thus it has zero 
mean, and a variance 



n(x,t)n(x',t')=a 2 6 n (x-x')6(t-t') . (99) 



(In a more general anisotropic situation, eq. (£JS|) has to be generalized to include different 
diffusivities v a along different directions. Such anisotropy is easily removed by rescaling 
the coordinates x a .) 

Eq. fl98|) can be regarded as the imaginary time Schrodinger equation for a particle 
in a random time dependent potential. It can be integrated to yield the continuous path 
integral 



W(x,t) = / Vx(t')exp - / dt'\--\ fi(x(t'),t' 

J (o,o) Jo V? Au 



(100) 
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describing the fluctuations of a directed polymer in a random medium (DPRM) [B31 . The 
71 th moment of W is computed by replicating the above path integral and averaging over 
n(x,t). It generalizes eq. ( [73p to 

W n ({x a },t)= [ ' Vx x tf)---Vx n {t') 



'({o},o) 
exp 




t ,„ -\Y, 6D ^{t')-x> P {t')) 



(ioi) 



with u oc a . The differential equation governing the evolution of W n (t) is, 

f~i 1/1/ T7 1 / 

-K^ = -TWn + ^J2 V ^ Wn+ 2^ 6D( - £a ~ Sf3)Wn -~ HnWn ■ (102) 

Evaluating the asymptotic behavior of W n (t) requires knowledge of the ground state 
energy of the Hamiltonian H n . Unfortunately, the exact dependence of the bound state 
energy on n is known only for D = (e oc n{n — 1)) and D = 1 (e oc n(n 2 — 1)). As 
discussed earlier, these two results can then be used to deduce the behavior of the bulk of 
the probability distribution for lnW(t). Elementary results from quantum mechanics tell 
us that an arbitrarily small attraction leads to the formation of a bound state in D < 2, 
but that a finite strength of the potential is needed to form a bound state in I? > 2. 
Thus, in the most interesting case of 2 + 1 dimensions we expect a non-trivial probability 
distribution, while the replica analysis provides no information on its behavior. In higher 
dimensions, there is a phase transition between weak and strong randomness regimes. For 



weak randomness there is no bound state and asymptotically W n (t) = W(t) n , indicating a 
sharp probability distribution. This statement has also been established by more rigorous 
methods [ p5| . There is another phase for strong randomness where the probability distri- 
bution for W(t) becomes broad. The resulting bound state has been analytically studied in 



a 1/D expansion valid for large D [26 1. The ground state wavefunction is rather complex, 
involving replica symmetry breaking. Note that the phase transition in the probability dis- 
tribution of the correlation function occurs in the high temperature phase of the random 
Ising model. The implications of this phase transition for bulk properties are not known. 
As the stiffness associated with line tension decreases on approaching the order/disorder 
phase transition of the Ising model, close to this transition the probability distribution for 
W(t) is likely to be broad. 
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As one of the simplest models of statistical mechanics in random systems (a "toy" 



spin glass), the problem of DPRM has generated considerable interest ||27|| . The model 
has been generalized to manifolds of arbitrary internal dimensions in random media |28|] , 
and treated by functional RG methods |29|]. The same model has also been studied by 



a variational approach that involves replica symmetry breaking [30|. The latter is also 
exact in the D — > oo limit. Directed paths have been examined on non-Euclidean lattices: 
In particular, the problem can be solved exactly on the Cayley tree ||31|| , where it has 
a transition between a "free" and a glassy state. There are also quite a few treatments 
based on a position space renormalization group scheme |H2] which becomes exact on a 
hierarchical lattice [ |55| . This lattice has no loops, and at the m + 1 th level is constructed 
by putting together 2 D branches, each containing two lattices of the m th level. Starting 
from a set of random bonds at the first level, the values of the sum W(m = log 2 t) are 
constructed recursively from 



W{m + 1,0) = ^2 W(m, ai)W(m, a 2 ) 



(103) 



a = \ 



where the greek indices are used to indicate specific bonds for a particular realization. 
Alternatively, these recursion relations can be used to study the evolution of the probability 
distribution for W [B4]]. The exponent uj m 0.30 for D = 1 is not too far off from the exact 
value of 1/3. 

Additional information about the higher dimensional DPRM is obtained by taking 
advantage of a mapping to the nonequilibrium problem of kinetic roughening of growing 
interfaces. Using the Cole-Hopf transformation |3J| , 



W{x, t) = exp 



Xh(x, t) 
2^ 



(104) 



eq. ([)8]) is transformed to the Kardar, Parisi, Zhang (KPZ) [|36| equation, 

M=^ + ^ 2 h--(Vh) 2 --ii(x,t) , (105) 

describing the fluctuations in height h(x, t) of a growing interface. A dynamical renormal- 
ization group (RG) analysis at the one-loop level J3"T| , |3"8| ] of this equation indicates that 
the effective coupling constant g = 4a 2 /u, satisfies the rescaling relation 



d JL = ( 2-D) 9+ K °V D - 3 '> g i 

d£ v )y D y 



(106) 
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where Kp is the D dimensional solid angle divided by (2n) D . The RG equation merely 
confirms the expectations based on the replica analysis: there is flow to strong coupling 
for D < 2, while there is a transition between weak and strong coupling behavior in higher 
dimensions. The RG equation has been recently extended to two loops |39|,(|0]]. According 



to one calculation [40|, there is no stable strong coupling fixed point for D > 2, and the 



nature of this phase remains a mystery. 



Since there are several comprehensive reviews of the KPZ equation [[4111 , I will not 
discuss its properties in any detail here. It suffices to say that there are many numerical 
models of growth that fall in the universality class of this equation. They are in complete 
agreement with the exactly known results for D = 1. The estimates for the exponent £ 



in higher dimensions are £ = 0.624 ± 0.001 for D = 2 g!J and C ~ 0.59 for D = 3 
The numerical results in higher dimensions are consistent with an exponent ( that gets 
closer to 1/2 as D — > oo. It is not presently known whether there is a finite upper critical 
dimension [p8| , |44}| beyond which £ = 1/2 exactly. 



I. Random Signs 

So far we focused on nearest neighbor bonds {Kij}, which though random, were all 
positive. For such couplings the ground state is uniform and ferromagnetic. The study 
of low temperature states is considerably more complicated for the random spin glass 
which describes a mixture of ferromagnetic and antiferromagnetic bonds. The competition 
between the bonds leads to frustration effects, resulting in quite complicated landscapes 
for the low energy states [Q. Here we shall explore the high temperature properties of spin 
glass models. To focus on the effects of the randomness in sign, we study a simple binary 
probability distribution in which negative and positive bonds of equal magnitude occur 
with probabilities p and 1 — p respectively. 

The computation of the high temperature series for the correlation function (along 
the diagonal) proceeds as before, and 

t 

W(x,t) = (a , o' :S ,t) = r t J2llvPi , (107) 

P i=l 

where r indicates the fixed magnitude of tanhX, while r\pi = ±1 are random signs. Since 
the elements of the sum can be both positive and negative, the first question is whether 
the system maintains a coherence in sign (at least for small p), i.e. what is the likelihood 
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that the two spins separated by a distance t have a preference to have the same sign. This 
question can be answered only in one and high dimensions. 

For the one dimensional chain the moments of W(t) are easily calculated as 

W^ = T nt x ni-2pY for n odd, (10g) 

I- 1 for n even. 

As all odd moments asymptotically decay to zero, at large distances W(t) is equally likely 
to be positive or negative. This is because the sign of the effective bond depends only 
on the product of the intermediate bonds and the possibility of a few negative bonds is 
sufficient to remove any information about the overall sign. From eq. ( p.08| ) , we can define 
a characteristic sign correlation length £ s = — l/ln(l — 2p). 



There is also a "mean-field" type of approach to the sign coherence problem ]45[ which 



is likely to be exact in high dimensions. For paths along the diagonal of the hyper cubic 
lattice, the mean value of W(t) is 



W(t) « [dr(l - 2p)] z . (109) 

Calculating the variance of W is complicated due to the previously encountered problem 
of intersecting paths. We can approximately evaluate it by considering a subset of paths 
contributing to the second moment as, 

P « [dr(l - 2p)] 2t + (olr 2 ) [dr(l - 2p)] 2(t_1) + (air 2 ) 2 [dr(l - 2p)] 2(t_2) + • • • + (dr 2 )* 

r2t [d(l-2p)f t+1) -d^ 

[d{l-2p)f-d 

(110) 

The first term in the above sum comes from two distinct paths between the end points; 

the second term from two paths that have their first step in common and then proceed 

independently. The m th term in the series describes two paths that take m steps together 

before becoming separated. The underlying assumption is that once the two paths have 

separated they will not come back together again. This Independent Path Approximation 

(IPA) is better justified in higher dimensions and leads to 

WHt) _ d{l-2pf-[d{l-2pfY t 

wit) 2 d{i-2pY-i ■ [ ] 

For small p such that d(l — 2p) 2 > 1, the above ratio converges to a constant as t — >• oo; 
the distribution is asymptotically sharp and the correlations preserve sign information. 
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However, if the concentration of negative bonds exceeds p c = ( 1 — l/yd) /2, the ratio 
diverges exponentially in t, indicating a broad distribution. This has been interpreted 
||46|| as signalling a "sign transition" . This argument suggests that there is a finite p c for 
all d > 1. However, it is important to note that the IPA ignores important correlations 
between the paths. Shapir and Wang [I7| criticize the assumption of independent paths and 
suggest that as intersections are important for d < 3, there should be no phase transition 
in these dimensions. However, the identification of the lower critical dimension for the sign 
transition is not completely settled. Numerical simulations based on the transfer matrix 
method for t of up to 600(11], as well as exact enumeration studies [IS| for t < 10, fail to 



find a phase transition in d = 2. The results suggest that if there is a phase transition in 
d = 2 it occurs for p c < 0.05. The phase diagram of a generalized model with complex 
phases has also been studied in higher dimensions ]^9| , |50[] . 

For p > p c the information on sign is lost beyond a coherence length £ s . If the system 
is coarse grained beyond this scale, the effective bonds are equally likely to be positive 
or negative. Thus we shall concentrate on the symmetric case of p = 1/2 in the rest of 



this section. This corresponds to the much studied ±J Ising spin glass [51] which will 



be discussed in more detail later on. We performed QT4 j transfer matrix computations on 



systems of up to size t = 2000, and averaged over 2000 realizations of randomness, on a 
VaxStationll. The random numbers (+1 or —1) were generated by a well tested random 



number generator [52 . Since W grows exponentially in t, In | W| has a well defined proba- 



bility distribution; we examined its mean In |W(£)|, and variance In |V^(t) | 2 — In |W(£)| 2 , 
for p — 1/2 (both signs equally probable). We also computed the typical excursions of the 
paths in the lateral direction as defined by 



■£ x xi\W(x,t)\' 



Hm - s x.\w<;,w ■ (112) 



and 



nr_ (Y.„4w{x,t)r- 



ix(t)] -[ %wm ) ' <n3) 

where [•]„„ denotes an average over the lateral coordinate at a fixed t, using a weight 

|^,i)| 2 . 

We confirmed that the average of In \W(t)\ is extensive (ln|W(t)| = (0.322 ± 0.001)*), 
while its fluctuations satisfy a power law growth t u , with uj = 0.33 ± 0.05. For several 
choices of t we also checked in detail that W(t) is positive or negative with equal probability. 
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For lateral excursions, we examined simulations with t = 4000, and with 200 realizations 
of randomness (reasonable data for fluctuations of In |W(£)| are only obtained from higher 



averaging). The results for [x 2 ] av and [x] 2 v appear to converge to a common asymptotic 
limit; fitted to a power law £ 2< > with £ = 0.68 ± 0.05. The scaling properties of |W(x,t)| 
thus appear identical to those of directed polymers with positive random weights! It should 
be noted, however, that using a similar procedure, Zhang [[53] concluded from fits to his 



numerical results a value of £ = 0.74 ± 0.01. Using a variety of theoretical arguments] 
he suggests u> = 1/2 and £ = 3/4. The exponent u> = 1/2 is clearly inconsistent with 



our data, while ( = 3/4 can be obtained if one fits only to [x] 2 v . Two subsequent, rather 



extensive, numerical studies [p4] , |55|1 shed more light on this problem. Both simulations 
seem to equivocally point to the importance of including corrections to scaling in the fits. 
In 1+1 dimensions they indeed find ui = 1/3 for the variance, and ( = 2/3 (with a large 
correction to scaling term) for transverse fluctuations. 

The similarity in the probability distributions of random weight and random sign 
problems can be understood by examination of the moments. The terms in W n correspond 
to the product of contributions from n independent paths. Upon averaging, if m paths 
cross a particular bond (0 < m < n), we obtain a factor of [1 + (— l) m ]/2, which is or 
1 depending on the parity of m. For odd n there must be bonds with m odd, and hence 



W 2n+1 = 0; which of course implies and follows from the symmetry p(W) = p(—W). 



For even moments W 2n , the only configurations that survive averaging are those in which 
the In replicated paths are arranged such that each bond is crossed an even number of 
times. The simplest configurations satisfying this constraint correspond to drawing n 
independent paths between the end points and assigning two replica indices to each. The 
above constraint is also satisfied by forming groups of four or higher even numbers, but such 
configurations are entropically unlikely and we shall henceforth consider only paired paths. 



There is an important subtlety in calculating W 2n from the n paired-paths: After two 
such paths cross, the outgoing pairs can either carry the same replica labels as the ingoing 
ones, or they can exchange one label (e.g. (12)(34) -» (12)(34), (13)(24), or (14)(23)). 
Therefore, after summing over all possible ways of labelling the paired paths, there is a 
multiplicity of three for each intersection. The n paired paths attract each other through 
the exchange of replica partners! 
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Although the origin of the attraction between paths is very different from the case of 
random weights, the final outcome is the same. The even moments in 1 + 1 dimension are 
related by an expression similar to eq. (|89|), 



lim W 2n {t) = W{t) 2 n exp [pn(n 2 - l)t] , (114) 

and the conclusions regarding In W(t) are the same as before. If, rather than having only 
one possible value for the magnitude of the random bond, we start with a symmetric 
distribution p(r), there will be an additional attraction between the paired paths coming 
from the variance of r 2 . This increases the bound state energy (and the factor p) in 
eq. ( |114| ) but does not affect the universal properties. 



J. Other Realizations of DPRM 

So far we focused on sums over DPRM as encountered in high temperature series 
of Ising models. In fact several other realizations of such paths have been discussed in 
the literature, and many more are likely to emerge in the future. One of the original 
motivations was to understand the domain wall of an Ising model in the presence of random 



bond impurities [jl|. As mentioned in the previous section, if all the random bonds are 
ferromagnetic, in the ground state all spins are up or down. Now consider a mixed state in 
which a domain wall is forced into the system by fixing spins at opposite edges of the lattice 
to + and — . Bonds are broken at the interface of the two domains, and the total energy of 
the defect is twice the sum of all the Ky crossed by the interface. In the Solid-On-Solid 
approximation, configurations of the domain wall are restricted to directed paths. The 
resulting partition function Z(t), can be computed by exactly the same transfer matrix 
method used to calculate W(t). Rather than looking at the finite temperature partition 
function, Huse and Henley [|T5| worked directly with the zero temperature configuration of 
the interface. 

Denoting by E(x, t) the minimum in the energy of all paths connecting (0, 0) to (x, t), 
oriented along the diagonal of the square lattice, it is possible to construct the recursion 
relation, 

E(x,t+l) = min{E(x-l,t)-2J x _ ht ,E(x + l,t)-2J x+1>t } , (115) 

closely related to eq. (|95"D . To find the actual configuration of the path, it is also necessary 
to store in memory one bit of information at each point (x,t), indicating whether the 

32 



minimum in eq. ( |115| ) comes from the first or second term. This bit of information indicates 
the direction of arrival for the optimal path at (x,t). After the recursion relations have been 
iterated forward to "time" step t, the optimal location is obtained from the minimum of the 
array {E(x, t)}. From this location the optimal path is reconstructed by stepping backward 



along the stored directions. This is how the pictures of optimal paths in refs. [|24] , |38H were 
constructed. These optimal paths have a beautiful ultrametric structure that resembles 
the deltas of river basins, and many other natural branching patterns. Finding the optimal 
interface is reminiscent of the travelling salesman problem. However, in this case, although 
the number of possible paths grow as 2*, their directed nature allows us to find the best 
solution in polynomial time. 

The statistics of the E(x, t) at T = are identical to those of In W(x, t): the optimal 
path wanders as t 2 / 3 , while the fluctuations in E(t) scale as t 1 / 3 [|T5| . It is frequently 
assumed that these fluctuations also set the scale of energy barriers that the interface must 
cross from one optimal state to another. Since such barriers grow with t, any activated 
process is slowed down to a logarithmic crawl |[15| . 

It has been suggested that optimal paths are also relevant to fracture and failure 



phenomena [56]. Imagine a two dimensional elastic medium with impurities, e.g. a network 



of springs of different strengths and extensions [57]. If the network is subjected to external 
shear, a complicated stress field is set up in the material. It is possible that non-linear 
effects in combination with randomness enhance the stress field along particular paths in 
the medium. Such bands of enhanced stress are visible by polarized light in a sheet of 
plexiglass. The localization of deformation is nicely demonstrated in a two dimensional 



packing of straws [ 58(1 - The roughness of the localization band is characterized by the 



exponent £ = 0.73±0.07, not inconsistent with the value of 2/3 for DPRM. The experiment 



was inspired by random fuse models |59| which apply a similar procedure to describe the 



failure of an electrical network. Hansen et al. |)(| suggest that at the threshold in all such 



models, failure occurs along an optimal path with statistics similar to a DPRM. Their 
numerical results obtain a roughness exponent of £ = 0.7 for the crack interface with a 
precision of about 10%. 

In fact, the minimal directed path was proposed in 1964 |H]] as a model for tensile 



rupture of paper. The variations in brightness of a piece of paper held in front of a light 
source are indicative of nonuniformities in local thickness and density p(x). Tydeman and 
Hiron suggested that rupture occurs along the weakest line for which the sum of p(x) is 
minimum. This is clearly just a continuum version of the optimal energy path in a random 
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medium. (Since the average of p(x) is positive, the optimal path will be directed.) This 
model was tested by Kertesz et al. |62| who used a tensile testing machine to gradually tear 



apart many sheets of paper. They found that the resulting rupture lines are self-affine, 
characterized by 0.63 < C < 0.72. 

The three dimensional DPRM was introduced |24| as a model for a polyelectrolyte in 
a gel matrix. Probably a better realization is provided by defect lines, such as dislocations 
or vortices, in a medium with impurities. There has been a resurgence of interest in this 
problem since it was realized that flux lines in high temperature ceramic superconductors 
are highly flexible, and easily pinned by the oxygen impurities that are usually present 
in these materials ||63| . Pinning by impurities is in fact crucial for any application, as 
otherwise the flux lines drift due to the Lorentz force giving rise to flux flow resistivity 



64 



K. Quantum Interference of Strongly Localized Electrons 

The wavefunctions for non-interacting electrons in a regular solid are extended Bloch 
states. In the presence of disorder and impurities, gradually more and more of these states 
become localized. This was first pointed out by Anderson |65| who studied a random 
tight-binding Hamiltonian 

H = J2 £ * a t a i + Yl V v a t a 3 ■ ( 116 ) 

Here e% are the site energies and Vij represent the nearest neighbor couplings or transfer 
terms. For simplicity we shall focus on 



v __ ( V if i, j are nearest neighbors 
lJ 1 otherwise 



otherwise 

so that all the randomness is in the site energies. This is just a discretized version of the 
continuum Hamiltonian H = z/V 2 + s{x), for a quantum particle in a random potential 
e(x). For a uniform e, the Hamiltonian is diagonalized by extended Fourier modes at = 
^2- exp (iq ■ x) a~/vN, resulting in a band of energies e(q) = e + 21/(cos gi + cos 02 + • ■ • + 
cos qa)- (The lattice spacing has been set to unity.) As long as the fermi energy falls within 
this band of excited states the system is metallic. 

In the random system the wave functions become distorted, and localized to the 
vicinity of low energy impurities ||65|| . This localization starts with the states at the edge 
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of the band and proceeds to include all states as randomness is increased. In fact in d < 2, 
as the diffusing path of a non-localized electron will always encounter an impurity, all states 



are localized by even weak randomness. The original ideas of Anderson localization[65|], and 
a heuristic scaling approach by Thouless |)6| , have been placed on more rigorous footing by 
perturbative RG studies [|67| -|69| . The perturbative approach emphasizes the importance of 
quantum interference effects in the weakly disordered metal. Weak localization phenomena 
include the effects of magnetic fields, spin-orbit (SO) scattering (corresponding respectively 
to interactions breaking time reversal and spin space symmetries) on the conductivity |7TJ , 
as well as predicting a universal value of the order of e 2 /h for conductance fluctuations 
|7l],[72[]. These phenomena can be traced to the quantum interference of time reversed 



paths in backs cattering loops and their suppression by magnetic fields and SO [fr3] : In 
the of absence SO, a magnetic field causes an increase in the localization length, and a 
factor of 2 decrease in the conductance fluctuations; with SO, it has the opposite effect of 



decreasing the localization length, while still reducing the conductance fluctuations [74,75 



An alternative description of these phenomena is based on the theory of random matrices 
79|, where the only input is the symmetries of the underlying Hamiltonian and their 



modification by a magnetic field. Mesoscopic devices at low temperature have provided 
many experimental verifications of weak localization theory |73,77| and there are many 
excellent reviews on the subject [[F0|,[F8 . 



When the electronic states at the fermi energy are localized, the material is an insulator 
and there is no conductivity at zero temperature. However, at finite temperatures there 
is a small conductivity that originates from the quantum tunneling of electrons between 



localized states, described by Mott's variable range hopping (VRH) process [ff9| : The 
probability for tunneling a distance t is the product of two factors 

p(t) ocexpf-— J xex p(-ir^) • ( 117 ) 

The first factor is the quantum tunneling probability and assumes that the overlap of the 
two localized states decays with a characteristic localization length £. The second factor 
recognizes that the different localized states must have different energies 5e (otherwise 
a new state is obtained by their mixture using degenerate perturbation theory). The 
difference in energy must be provided by inelastic processes such as phonon scattering, 
and is governed by the Boltzmann weight at temperature T. The most likely tunneling 
sites must be close in energy. If there is a uniform density of states N(ef) in the vicinity 
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of the fermi energy, there are roughly N(sf)t d candidate states in a volume of linear size 
t in d dimensions, with the smallest energy difference of the order of 5e oc [N(ef)t d ) 
Thus the two exponential factors in eq. (|117| ) oppose each other, encouraging the electron 
to travel shorter and longer distances respectively. The optimal distance scales as 

t^(T /T)^ , (118) 

with To oc (kBN(ef)t; d ) , diverging at zero temperature. 

In the strongly localized regime, the optimal hopping length is many times greater 
than the localization length £. The localized sites are then assumed to be connected by a 



classical random resistor network |8JJ . Since the individual resistors are taken from a very 



wide distribution, it is then argued pi] that the resistance of the whole sample is governed 



by the critical resistor that makes the network percolate. This leads to a dependence 

a(T)=a exp[-(T /T)^T] , (119) 

for the conductivity. This behavior has been verified experimentally both in two and 
three dimensions |82|J. Due to the difficulty of measuring variations in the much smaller 



conductivities of insulators, there have been relatively few studies of the conductivity and 



its fluctuations for strongly localized electrons. Nonetheless, recent experiments [S3] find a 
positive MC in S'z-inversion layers, GaAs and In^O^-x films. Furthermore, the observed 
reproducible conductance fluctuations are quite suggestive of quantum interference (QI) 
effects. However, the magnitudes of these fluctuations grow with lowering temperature, 
and are about 100 times larger than e 2 /h at the lowest temperature. 

Clearly a different theory is needed to account for QI effects in the strong localization 
regime. The most natural candidate is the quantum overlap factor in eq. (|117|) . Nguyen, 
Spivak, and Shklovskii (NSS) have proposed a model that accounts for QI of multiply 
scattered tunneling paths in the hopping probability: In between the phonon assisted 
tunneling events the electron preserves its phase memory. However, at low temperatures it 
tunnels over very large distances according to eq. ( |118|) , and encounters many impurities. 



The overall tunneling amplitude is then obtained from the sum over all trajectories between 
the initial and final sites. NSS emphasized that since the contribution of each trajectory 
is exponentially small in its length, the dominant contributions to the sum come from 
the shortest or forward scattering paths. The traditional explanations of weak localization 
phenomena which rely on the QI of back scattering paths are therefore inappropriate to 
this regime. This picture is clearly reminiscent of the directed paths and will be developed 
more formally in the next section. 
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L. The Locator Expansion and Forward Scattering Paths 

The overlaps in the insulating regime can be studied by performing a "locator" expan- 
sion ||65|| ; valid in the limit \V%j\ = V <C {E — Ei), where E is the electron energy. Indeed, for 
V = 0, the eigenf unctions are just the single site states, and the localization length is zero 
(no transfer term). For V/(E — Si) <C 1, various quantities can be obtained perturbatively 
around this solution, as expressed by the Lippman-Schwinger equation |8|[ 



l* + ) = l$)+ p I M -^ + ) ■ (120) 

E — Hq + io 



The bare Hamiltonian 



• + a 7 - 



H = J2 £ i a i 

i 

has no nearest-neighbor coupling, while the perturbation 

V = Yl Vi J a t a J > 

describes the small transfer terms. |$) represents the state with a localized electron at the 
initial site (or incident wave), |^/ + ) the state where a localized electron is at the final site. 
In the coordinate representation, the wavefunctions are exponentially localized around the 
impurity sites and there are no propagating waves since electrons can only tunnel under a 
potential barrier. (This situation was first addressed in detail by Lifshits and Kirpichenko 
|85|| .) We can now iterate this implicit equation to obtain an expansion in powers of the 
ratio V / (E — Ei) as 

|W+) = |$) + V|$) + V V|$) + --- . (121) 

1 ' ' ' E-H + i5 ' ' E-H +iS E-H + i5 ' ' v ; 

Acting with (\I/ + | on the left and taking S to zero, we obtain the overlap between the two 
states 

^+|$+) = ^+|$) + ^+|^_V|*> + <*+|_J_V-J— V|*) + --- . (122) 

Hi — £±0 Hi — £±0 Hi — 11Q 

For a more general transfer term V connecting all sites, the first term represents an electron 
starting from the initial site and ending at the final site without scattering (the overlap 
(\I/ + |$)); the second term represents electrons scattering once off intermediate sites, the 
third, scattering twice, etc.. The operator V acting on |<&) produces a factor V for each 
segment crossed, and Hq acting on a particular site i results in £;, the bare site energy. 
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Thus we finally arrive at a simple expression for the amplitude or the Green's function 
between the initial and final states as 

(9+\9+) = ($\G(E)\* + ) = VJ2I[ 1 ^—- ■ (123) 

The terms in the above perturbation series correspond to all paths F connecting the 
end points; i-p label the sites along each path. Except that the random variables appear 
on the sites rather than the bonds of the lattice, this sum over paths is quite reminiscent 
of the corresponding one for the correlation functions of the random bond Ising model. 
There is, however, one complication that distinguishes the localization problem: The en- 
ergy denominators in eq. ( |123|) may accidentally be zero, invalidating the perturbation 
series. Physically, this corresponds to intermediate sites that are at the same energy as the 
external points. Presumably in this case a degenerate perturbation theory has to be used 



to construct the wavefunction. NSS |46|] circumvent this issue by considering initial and 
final sites of approximately the same energy ep = E = 0, while the intermediate sites have 
energies £{ = ±U with equal probability. All the energy denominators in eq. ( |123| ) now 
contribute the same finite magnitude U, but random signs n ir = e ir /U . The justification is 
that the Mott argument implicitly assumes that the lowest energy 5e occurs at a distance 
t, and that there are no intermediate sites that are more favorable. However, it is not 
clear that due to the very same considerations, we should not include some dependence of 
the effective energy gap U on t. We shall set aside such considerations and focus on the 
properties of the NSS model in the remainder. 

A path of length £ now contributes an amplitude U(V/U) to the sum, as well as an 
overall sign. In the localized regime the sum is rapidly convergent, dominated by its lowest 



order terms |65| . In general, the sum is bounded by one in which all terms make a positive 
contribution, i.e. by a lattice random walk which is convergent for z(V/U) < 1, where 
z is the lattice coordination number. This provides a lower bound for the derealization 
transition, and the series is certainly convergent for smaller values of V/U. As in the 
Ising model we expect loops to become important only after the transition, while in the 
localized phase typical paths are directed beyond the localization length £. For (V/U) <C 1, 
the localization length is less than a single lattice spacing, and only directed (forward 
scattering) paths need to be considered. Loops ( back scattering paths) are irrelevant in the 
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renormalization group sense. For sites separated by a distance t + 1 along a diagonal of 
the square lattice, eq. ( |123|) is now simplified to 



' P j=l 



(124) 



which is identical to eq. (|f07| ) with (V/U) replacing r. The diagonal geometry maximizes 
possible interference by having a large number of shortest paths. For tunneling along the 
axes rather than the diagonal of a square lattice there is only one shortest path. Then, 
including longer paths with kinks is essential to the interference phenomena. However, 
the analogy to previous results suggests that the universal behavior is the same in the two 
cases while the approach to asymptotic behavior is much slower in the latter. 
Using the equivalence to eq. ( |114j ), in conjunction with eq. (|93"D, results in 



lim In |<z|C|/>| 2 

t — >oo 



In 



< 



t _ /Ot = _2t(e - 1 +C 1 ) 



(125) 



where we have defined local and global contributions to the effective localization length, 
respectively given by 



Co 



-a) 



-i -i 



and f _ 1 



P 



(126) 



The QI information is encoded in 2£ l = p. Numerical estimates indicate that for the 
NSS model £ g ~ 40, and confirm that the width of the distribution scales as 



<51n|(W>l 



t 



1/3 



(127) 



Since t oc T -1 / 3 in Mott VRH, we expect fluctuations in log-conductivity to grow as 
T -1 ' 9 for T — >• 0, in qualitative agreement with the experimental results of ref. ||83|| . (A 
quantitative test of this dependence has not yet been performed.) 



M. Magnetic Field Response 

All that is needed to include a magnetic field B in the tight binding Hamiltonian of 
eq. ( |116| ) is to multiply the transfer elements Vij by exp (Aij), where A^ is the line integral 
of the gauge field along the bond from i to j. Due to these factors, the Hamiltonian becomes 
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complex and is no longer time reversal symmetric (H* ^ H). In the parlance of random 
matrix theory |76j , the Hamiltonian with B = belongs to the orthogonal ensemble, while 



a finite field places it in the unitary ensemble. Actually, random matrix theory recognizes 
a third (symplectic) ensemble of Hamiltonians which are time reversal symmetric, but not 
invariant under rotations in spin space. Up to this point we had not mentioned the spin 
of the electron: The states of eq. (|116| ) are thus doubly degenerate and can be occupied 
by (non-interacting) up or down spin states. We can remove this degeneracy by including 
spin-orbit (SO) scattering, which rotates the spin of the electron as it moves through the 
lattice. 

The generalized tight binding Hamiltonian that includes both the effects of SO scat- 
tering and magnetic field is 

H = J2 £ * a L a ^ + J2 v iJ,™' eiAii a>t* a 3,<r' ■ ( 128 ) 

i,cr <ij>,crcr / 

The constant, nearest- neighbor only hopping, elements V in eq. ( |116| ) are no longer diag- 
onal in spin space. Instead, each is multiplied by Uij, a randomly chosen 577(2) matrix 
which describes the spin rotation due to strong SO scatterers on each bond [f76| . Eq. (|123| ) 
for the overlap of wavefunctions at the two end-points must now include the initial and 
final spins, and eq. ( |124|) for the sum of directed paths generalizes to 



t 
A = (i*\G(0)\f*') = V(V/U)*J(t) ; J(t) =YsH VP 3 e lAp ^^Up h P {3+ i) • (129) 

P 3=1 

After averaging over the initial spin, and summing over the final spin, the tunneling prob- 
ability is 

T=±tr(AU) = V\V/U) 2t I(t) ; I(t) = \ tr(J+J) . (130) 

We numerically studied the statistical properties of I(t), using a transfer matrix 
method to exactly calculate I up to t = 1000, for over 2000 realizations of the random 
Hamiltonian. We again found that the distribution is broad (almost log-normal), and 
that the appropriate variable to consider is ln/(t). In all cases the mean of ln/(t) scaled 
linearly with t, while its fluctuations scaled as t w with uj w 1/3 [jS6H38f. For the sake of 
comparison with experiments we define a log-magneto conductance (MC) by 



MC{t, B) = lnl(t, B)- In J(t, 0) . (131) 
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We find numerically that the magnetic field always causes an enhancement in tunneling (a 
positive MC), but that the asymptotic behavior is quite distinct in the presence or absence 
of SO scattering. 

(1) In the absence of SO, the MC is unbounded and grows linearly with t. This can be 
interpreted as an increase in the global contribution to the localization length. The 
numerical results indicate that for small B, the change in slope is proportional to 
S 1 / 2 . Indeed the data for different t and B can be collapsed together, using the fit 

/ , \ 1/2 

MC(t,B) = (0.15 ±0.03) (J- J t , (132) 

where <p = Ba 2 is the flux per plaquette, and <fio is the elementary flux quantum. 

(2) In the presence of SO, the MC quickly saturates with t and there is no change in the 
localization length. The data can still be collapsed, but by using Bt 3 / 2 as the scaling 
argument, and we find 

MC SO (t,B) = {f/ o25 »*£<! . (133) 

We can gain some analytic understanding of the distribution function for I(t, B) by 



examining the moments I(t) n . From eqs. ( |129| ) and ( |130|) we see that each I(t) represents 



a forward path from i to /, and a time reversed path from / to i. For I(t) n , we have to 
average over the contributions of n such pairs of paths. Averaging over the random signs of 
the site energies forces a pairing of the In paths (since any site crossed by an odd number of 
paths leads to a zero contribution) [J5BJ. To understand the MC, it is useful to distinguish 
two classes of pairings: (1) Neutral paths in which one member is selected from J and the 
other from jt. Such pairs do not feel the field since the phase factors of e lA picked up 
by one member on each bond are canceled by the conjugate factors e~ lA collected by its 
partner. (2) Charged paths in which both elements are taken from J or from J'. Such 
pairs couple to the magnetic field like particles of charge ±2e. In the presence of SO, we 
must also average over the random SU(2) matrices. From the orthogonality relation for 
group representations |89|1 , we have 

/AAA r 

r fc (fir)*,r fc (g)i'j'W(ai, ■■■, a n )dai ■ ■ ■ da n = ^ / W(ati, • • • , a n )dai ■ ■ ■ da n , 

(134) 

where T k (g)ij is the ij matrix element of a representation of the group element g, 

W(ai, ■ • • ,a n ) is an appropriate weight function so that the matrix space is sampled 
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uniformly as the continuous parameters a±, ■ ■ ■ ,a n vary (e.g. Euler angles for a represen- 
tation of SU(2)). Finally A& is the order of the representation k. Choosing the Euler angle 
parametrization of SU{2) it can be shown that the only nonzero paired averages are 



1 .. .. 1 .. .. 1 



UapW aP = - , W T tW U = 2 » U n U^ = -- , (135) 

and their complex conjugates. Thus SO averaging forces neutral paths to carry parallel 
spins, while the spins on the two partners of charged paths must be antiparallel. 

We next consider the statistical weights associated with the intersections of paths. 
These weights depend crucially on the symmetries of the Hamiltonian in eq. ( |128| ): For 
B = and without SO, the Hamiltonian has orthogonal symmetry. All pairings are allowed 
and the attraction factor is 3, since an incoming (12)(34) can go to (12)(34), (13)(24), 
or (14) (23). Note that even if both incoming paths are neutral, one of the exchanged 
configurations is charged. A magnetic field breaks time reversal symmetry, discourages 
charged configurations, and reduces the exchange attraction. The limiting case of a 'large' 
magnetic field is mimicked by replacing the gauge factors with random phases. In this 
extreme, the Hamiltonian has unitary symmetry and only neutral paths are allowed. The 
exchange factor is now reduced to 2; from (11*)(22*) -> (11*)(22*), or (12*)(21*). 
With SO averaging, we must also take into account the allowed spin exchanges: Two 
neutral paths entering the intersection can have indices (act), (eta) or (act), (act); there are 
2 possibilities for the first (a =] or |) and two for the second (a) is antiparallel to a). In 
the former case, however, there are two exchanges preserving neutrality, while in the latter 
only one exchange is possible satisfying this constraint. Hence an overall multiplicity of 
[2 x 2 + 2 x 1] x (1/2) 2 = 3/2 is obtained, where the (1/2) 2 comes from the averages in 
eq. ( |135| ) . Thus the intersection of two paired paths results in an exchange attraction of 
3/2; a signature of the symplectic symmetry. 

Based on the above symmetry dependent statistical attraction factors, we can provide 
an understanding of the numerical results for MO The sum over n attracting paths again 
leads to 

(I(t) n ) = A(n)2 nt exp[pn(n 2 - l)t] , (136) 

where we have also included an overall n-dependent amplitude. Without SO, the magnetic 
field gradually reduces the attraction factor from 3 to 2 leading to the increase in slope. 
Addition of SO to the Hamiltonian has the effect of suddenly decreasing the attraction 
factor to 3/2. Why does the addition of the magnetic field lead to no further change in p in 
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the presence of SO? Without SO, the origin of the continuous change in the attraction factor 
is a charged bubble that may appear in between successive intersections of two neutral 
paths. In the presence of SO, from the averages in eq. ( |135| ) we find the contribution of such 



configurations to be zero. To produce intermediate charged paths (with their antiparallel 
spins), the entering pair must have indices of the type (ii), (ii) (where j =|, and f =!)• 
Within the bubble we can have intermediate sites labeled (jj) and (kk) which must be 
summed over due to matrix contractions. It is easy to check that, independent of the choice 
of j, if the incoming and outgoing spins (i and to) are the same on a branch it contributes 
a positive sign, while if they are opposite the overall sign is negative. However, for any 
choice of i and to, one may choose similar (e.g. i — > m on both branches), or opposite (e.g. 
i —>■ m on top and i-»mon lower branch) connections. The difference in sign between the 
two choices thus cancels their overall contributions. Hence the neutral paths traverse the 
system without being affected by charged segments. In a magnetic field, their attraction 
factor stays at 3/2 and p = t;^ 1 is unchanged. 

The smaller positive MC observed in the simulations is due to changes in the amplitude 
A{n) in eq. (|136|) . This originates from the charged paths that contribute to tunneling at 
small B but are quenched at higher B. However, due to their lack of interactions, we may 
treat the charged and neutral paths as independent. At zero field any of the pairings into 
charged and neutral paths is acceptable, while at finite fields only neutral pairs survive. 
This leads to a reduction in the amplitude A(n) for n > 2, but an increase in In I (a 
positive MC). The typical value of In I thus increases by a t independent amount. This 
behavior is similar to the predictions of IPA, and is indeed due to the independence of 
charged and neutral paths. Since the typical scale of decay for charged paths depends on 
the combination Bt 3 ' 2 (typical flux through a random walk of length t) , we can explain 
the scaling obtained numerically in eq. ( |133|) . 

The exchange attraction between neutral paths can also be computed for (unphysical) 
SU(n) impurities and equals 1 + 1/n, which reproduces 2 for U(l) or random phases, and 
3/2 for SU(2) or SO scattering. The attraction vanishes in the n — > oo limit, where the 
paths become independent. The statistical exchange factors are thus universal numbers, 
simply related to the symmetries of the underlying Hamiltonian. The attractions in turn 
are responsible for the formation of bound states in replica space, and the universal scaling 
of the moments in eq. ( |136| ). In fact, since the single parameter p completely characterizes 
the distribution, the variations in the mean and variance of ln/(t) should be perfectly 
correlated. This can be tested numerically by examining respectively coefficients of the 
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mean and the variance for different cases. All results do indeed fall on a single line, 
parametrized by p. The largest value corresponds to the NSS model for B = and no SO 
(orthogonal symmetry, exchange attraction 3). Introduction of a field gradually reduces p 
until saturated at the limit of random phases (unitary symmetry, exchange attraction 2). 
SO scattering reduces p further (symplectic symmetry, exchange attraction 3/2). 

N. Unitary Propagation 

We can put together the results discussed so far by generalizing eq. ( |98| ) to allow for 
complex (and matrix valued) parameters. In the originally encountered directed polymer 
(DP), the parameters v > and p appearing in this equation were both real. To discuss 
the wavefunction in a magnetic field, we have to allow p to take complex values. Finally, 
SO scattering is included by generalizing W to a two component spinor, and using matrix 
valued p. We found that in all these cases the statistical behavior of In VF(x, t) is the same. 
Is this a general property of eq. (|98]), independent of the choice of parameters? A special 
limit of this equation is when both p — > —ip and v — > — iv are purely imaginary. Then 
eq. (^D reduces to the Schrodinger equation 

t^=[vV 2 +p(x,t)]W , (137) 

for a particle in a random time dependent potential. This equation has been considered in 
the context of particle diffusion in crystals at finite temperature p0^p2|, and to model the 
environment of a light test particle in a gas of much heavier particles |)3]. Several authors 



9^H96[1 have also suggested that the diffusion of directed wave fronts in disordered media 



are described by eq. ( |137| ). 

The path-integral solution to eq. ( |137| ) is 



w <* t] = C vAt,) exp f -' SI * fe 80 2 + Mx(t ' * 



(138) 



where x(t') now describes a path in d— 1 dimensions. In writing eq. ( |138|) , we have chosen 
the standard initial condition that at time t = 0, the "wave function" is localized at the 



origin. The beam positions <x 2 > and <x> 2 characterize the transverse fluctuations of 
a directed beam W about the forward path of least scattering. Here we use < • • ■ > to 
indicate an average with the weight |VF(a;,£)| 2 for a given realization, and 7TT to indicate 



quenched averaging over all realizations of randomness. Roughly speaking, <x> 2 describes 
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the wandering of the beam center, while <x 2 > — <x> 2 provides a measure of the beam 
width. 

A special property of eq. (|137|) which is valid only for real v and \i is unitarity, i.e. 
the norm J <ix| VF(x, t)\ 2 is preserved at all times. (In the DP and tunnelling problems, the 
norm clearly decays with the length t.) This additional conservation law sets apart the 
random directed wave problem from DP, and in a sense makes its solution more tractable. 
This unitarity is of course a natural consequence of particle conservation for the Schrodinger 
equation, but has no counterpart for directed wave propagation. It is likely that a beam 
of light propagating in a random medium will suffer a loss of intensity, due to either 
back-reflection, inelastic scattering, or localization phenomena |)8 . 



A number of efforts at understanding unitary propagation in random media have 



focused on the scaling behavior of the beam positions <x 2 > and <x> 2 at large t. Lattice 
models have been used here with some success. It has been shown using density-matrix 



techniques, for instance, that <x 2 > scales linearly in time as a consequence of unitarity [^TJ: 



recent numerical simulations |]99| , |100| also support this view. The scaling behavior of <x> 2 
at large t, however, is somewhat more complicated. An early numerical simulation in this 
area in ref. |)6| , employed a discretization procedure in which the norm of the wave function 
was not strictly preserved. In 2d, they found that |<x>| grew superdiffusively as £*> with 
C ~ 3/4, while in 3d, they found a phase transition separating regimes of weak and strong 
disorder. Subsequent numerical studies |)I| on directed waves in 2d cast doubt on the 



validity of these results when the time evolution is strictly unitary, indicating that <x> 2 
scales subdiffusively in 2d with C ps 0.3. 

Somewhat surprising is the fact that a continuum formulation of the wave problem 
leads to different results. An exact treatment of the continuum Schrodinger equation(|137Q 



has been given by Jayannavar and Kumar |j92|| . They show that for a random potential 
5-correlated in time, <x 2 > ~ t 3 as t — > oo. This behavior is modified when there are short- 



range correlations in time|93|], but the motion remains non-diffusive in that the particle is 
accelerated indefinitely as t — > oo. Lattice models introduce a momentum cutoff p m ax ~ 
a -1 , where a is the lattice spacing, and therefore do not exhibit this effect. The momentum 
cutoff generated by the lattice discretization is in some sense artificial. Nevertheless, in a 
real fluctuating medium, we do expect on large time scales to recover the lattice result, i.e. 
normal diffusion. The reason is that dissipative effects do generate an effective momentum 
cutoff in most physical systems. (Strictly speaking, even in the absence of dissipation, 
relativistic constraints lead to a velocity cutoff v = c.) The presence of such a cutoff for 
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the wave propagation problem, and hence the physical relevance of lattice versus continuum 
models, is still a matter of debate. While there is no underlying lattice, one suspects on 
physical grounds that there does exist an effective momentum cutoff for propagating waves, 
related to the speed of light in the background medium. 

Previous numerical investigations of this problem start with a discretization of the 
parabolic wave equation in eq. ( |137|) . By contrast, we assume that the path integral 



representation is more fundamental and provide a direct discretization of eq. ( |138|) that 



preserves unitarity [[101| . For concreteness, we introduce the model in 2d. A discussion of 



its generalization to higher dimensions is taken up later. As usual, we identify the time axis 
with the primary direction of propagation and orient it along the diagonal of the square 
lattice. The wave function is defined on the bonds of this lattice. We use W±(x, i) to refer 
to the amplitude for arriving at the site (x,t) from the ±x direction. At t = 0, the wave 
function is localized at the origin, with W±(0, 0) = l/v2. Transfer matrix techniques are 
now used to simulate diffusion in the presence of disorder. At time t, we imagine that 
a random scattering event occurs at each site on the lattice at which either W+(x,t) or 
W- (x, t) is non-zero. We implement these events by assigning to each scattering site a 2 x 2 
unitary matrix S(x,t). The values of the wave function at time t + 1 are then computed 
from the recursion relation: 

W-(x + l,t+l)\ _ fSn(x,t) S 12 (x,t)\ fW-(x,t) 



W+(x-l,t+l)J \S 21 (x,t) S 22 (x,t)J \W+{x,t)J (139) 

The S'-matrices are required to be unitary in order to locally preserve the norm of the wave 
function. As a particular realization, we may consider the rotation matrix 

o (f) ,«-( cos ^/ 2 ) e%4, sin ^/ 2 ) e ~ i4> \ rum 

b ^^)-{-sin(6/2)e^ cos (9/2) e~* * ) ' {14U) 

A physical realization of this model is obtained by placing semi polished mirrors of variable 
thickness, parallel to the t axis, on the sites of a square lattice. Within this framework, it 
should be clear that the value of W± (x, t) is obtained by summing the individual amplitudes 
of all directed paths which start at the origin and arrive at the point (x, t) from the ±x 
direction. We thus have a unitary discretization of the path integral in eq. (|138| ) in which 
the phase change from the potential n(x, t) is replaced by an element of the matrix S(x, t). 
A lattice S'-matrix approach for the study of electron localization and the quantum Hall 



effect has been used by Chalker and Coddington ||102|| . A related model has also been 



recently proposed ||103|| to investigate the localization of wave packets in random media. 



These models also include back scattering and hence involve a larger matrix at each site. 
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O. Unitary Averages 

A particularly nice feature of unitary propagation is that the weights W(x,t) are 
automatically normalized. In particular, we are interested in the beam positions 

<x 2 (t)> = J2 P ( x ^) x 2 > ( 141 ) 

x 

and 

<x(t)> 2 = Y^ p ( x i^) p ( x 2,t) xix 2 , (142) 

X\,X2 

where P(x,t) is the probability distribution function on the lattice at time t, defined by 

P(x,t)=\ W+(x,t) | 2 + | W-(x,t) | 2 . (143) 

(Defining the weights directly on the bonds does not substantially change the results.) 
Note that unlike the DP problem, P(x,t) is properly normalized, i.e. 

Y,P(x,t) = l , 



and eqs. (|141|) and ( |142|) are not divided by normalizations such as ^ x P(x,t). This 

simplification is a consequence of unitarity and makes the directed wave problem tractable. 

The average 7TT , in eqs. ( 141 ) and ( 142|) is to be performed over a distribution of 



S'-matrices that closely resembles the corresponding distribution for \i in the continuum 
problem. However, by analogy to the DP problem we expect any disorder to be relevant. 
Hence, to obtain the asymptotic scaling behavior, we consider the extreme limit of strong 
scattering in which each matrix S(x,t) is an independently chosen random element of the 
group U(2). With such a distribution we lose any pre-asymptotic behavior associated 



with weak scattering |J3|. The results are expected to be valid over a range of length scales 
a<x<^, where a is a distance over which the change of phase due to randomness is 
around 2-7T, and £ is the characteristic length for the decay of intensity and breakdown of 



unitarity. In the language of path integrals, the quantity P(x, t) represents the average 
over a conjugate pair of paths (from W± and W± respectively.) As in the random sign 



problem, the paths must be exactly paired to make a non-zero contribution (since S aj 3 = 0. 
In the strong disorder limit, each step along the paired paths contributes a factor of 1/2. 
(It can be easily checked from eq. (PO) ) that \S af3 \ 2 = cos 2 {9/2) = sin 2 (#/2) = 1/2.) Thus, 
in this limit, the effect of an impurity at (x, t) is to redistribute the incident probability 
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flux P(x,t) at random in the +x and — x directions. On average, the flux is scattered 

symmetrically so that the disorder-averaged probability describes the event space of a 

classical random walk, i.e. 

t\ 



P(x,t) = —. — ",, N , . (144) 

( t—x \ i f t+x \ i v / 



Substituting this into eq. ( |141| ), we find <x 2 (t)> = t, in agreement with previous stud- 
Hi- 



les 



Consider now the position of the beam center <x(t)> 2 , given by eq. (|142| ) . Unlike 



P(x,t), the correlation function P(xi,t)P(x 2 ,t) does not have a simple form. It involves 
a sum over four paths, collapsed into two pairs by randomness averaging. The center of 
mass coordinate R = (x% + x 2 )/2, performs a random walk with R 2 = t/2. Let us define 
a new correlation function for the relative coordinate r = x 2 — Xx, as 

W 2 (r, t) = J2 P(R ~ r/2, t)P(R + r/2, t) , (145) 

R 

with the initial condition 

W 2 (r : t = 0)=5 rfi . (146) 

The value of W2(r, t) is the disorder-averaged probability that two paired paths, evolved 
in the same realization of randomness, are separated by a distance r at time t, and can 
be computed as a sum over all configurations that meet this criteria. Consider now the 
evolution of two such pairs from time t to t + 1. Clearly, at times when r /0, the two 
pairs behave as independent random walks. On the other hand, when r = 0, there is an 
increased probability that the paths move together as a result of participating in the same 
scattering event. An event in which the pairs stay together is enhanced (since \S a /3\ = 



cos 4 (6/2) = sin (9/2) = 3/8), while one in which the pairs separate is diminished (since 
sin (0/2) cos 2 (9/2) = 1/8). These observations lead to the following recursion relation for 
the evolution of W 2 (r, t) , 

W 2 (r,t + 1) = (i±^) W ,{r,t)+(^*\ W 3 {r-2,t)+(?-^=*) W 2 (r + 2,t). 



(147) 
The parameter e > measures the tendency of the paths to stick together on contact. 
(If the S'-matrix is uniformly distributed over the group U(2), then e = 1/4.) Note that 
J2 r W 2 (r) is preserved, as required by unitarity. 
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Using eq. ( |147| ) , we evolved W% (r, t) numerically for various values of < e < 1 up 
to t < 15000. The position of the beam center was then calculated from 

<x(t)>2 = J2 (r 2 - j) P(R-r/2,t)P(R + r/2,t) = l - - 1^ W 2 (r,t) r 2 . (148) 



The results suggest quite unambiguously that <x(t)> 2 scales as £ 2< =, with £ = 1/4. We 
emphasize here the utility of the S'-matrix model for directed waves in random media. 
Not only does our final algorithm compute averages over disorder in an exact way, but 
it requires substantially less time to do so than simulations which perform averages by 
statistical sampling as in DPRM. We have in fact confirmed our 2d results with these 
slower methods on smaller lattices (t < 2000). 

The model is easily extended to higher dimensions. The wave function takes its values 
on the bonds of a lattice in d dimensions. Random d x d dimensional S'-matrices are then 
used to simulate scattering events at the sites of the lattice. When the matrices 5'(x, t) are 
distributed uniformly over the group U(d), the same considerations as before permit one 
to perform averages over disorder in an exact way. In addition, one obtains the general 



result for d > 2 that <x 2 > scales linearly in time. The computation of <x> 2 in d > 2, 



of course, requires significantly more computer resources. We have computed <x> 2 on a 
d = 3 body-centered cubic lattice, starting from the appropriate generalization of eq. ( |147|) . 



The results for t < 3000, indicate that <x> 2 scales logarithmically in time. 

The above numerical results can be understood by appealing to some well-known 
properties of random walks. Consider a random walker on a D = d — 1 dimensional 
hypercubic lattice. We suppose, as usual, that the walker starts out at the origin, and that 
at times t = 0, 1, 2, • • •, the walker has probability < p < 1/2D to move one step in any 
lattice direction and probability 1 — 2Dp to pause for a rest. The mean time to spent by 
the walker at the origin grows as ||101|| 



£2 (-0 = 1) 

/,-. - { hit (-0 = 2) . (149) 

constant {D = 3) 



The numerical results indicate a similar scaling for the wandering of the beam center <x> 2 
in d = D + 1 dimensions, for d = 2 and d = 3. We now show that this equivalence is not 
coincidental; moreover, it strongly suggests that d u = 3 is a critical dimension for directed 
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waves in random media. To this end, let us consider a continuum version of eq. ( |147|) , 
which in general dimensions takes the form 

W 2 (r,t + l)=W 2 (r,t) + V 2 [W 2 (l-e5 D (r) + ■••)] . (150) 

The asymptotic solution for e = is just a gaussian packet of width r 2 = 2t. We can next 
perform a perturbative calculation in e. However, simple dimensional analysis shows the 
corrections scale as powers of e/r D ~ et~ D ' 2 , and thus 

f 1 ^^( r -*) = (S^-p(-s)[ 1 + (^ ( "- 1,/2 )] ■ ( 1B1 > 

Applying the above results to the continuum version of eq. (|147|) , gives 

r 
-!_!/■ d D rV 2 [^ ^ _ e6 D {r) + . . .)] (152) 

~ l --\J d D vW 2 (l-e5 D (v))=eW 2 {0 1 t) . 
Summing both sides of this equation over t, one finds 

<x(t + l)> 2 = eVH/ 2 (0,t , )~ / dt'(4nt')- D / 2 . 



(153) 



The final integral is proportional to the time a random walker spends at the origin, and 
reproduces the results in eq. ( |149| ) . 

We can also regard W 2 (r, t) as a probability distribution function for the relative 
coordinate between two interacting random walkers. In this interpretation, the value of 
e in eq. ( |147| ) parametrizes the strength of a contact interaction between the walkers. If 
e = 0, the walkers do not interact at all; if e = 1, the walkers bind on contact. According to 



eq. ( |153| ), the wandering of the beam center <x(t)> 2 is proportional to the mean number 



of times that the paths of these walkers intersect during time t. If e = 0, the number of 
intersections during time t obeys the scaling law in eq. ( |149| ), since in this case, the relative 



coordinate between the walkers performs a simple random walk. Numerical results indicate 
that the same scaling law applies when < e < 1: the contact attraction does not affect 
the asymptotic properties of the random walk. In summary, three classes of behavior are 
possible in this model. For e = 0, i.e. no randomness, the incoming beam stays centered 
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at the origin, while its width grows diffusively. For < e < 1, the beam center, <x> 2 , 
also fluctuates, but with a dimension dependent behavior as in eq. ( |149|) . In the limit of 



e = 1, interference phenomena disappear completely. In this case, the beam width is zero 
and the beam center performs a simple random walk. 

We conclude by comparing the situation here to that of DPRM. In the replica approach 
to DPRM, the n-th moment of the weight W(x,t) is obtained from the statistics of n 
directed paths. Disorder-averaging produces an attractive interaction between these paths 
with the result that they may form a bound state. In d < 2, any amount of randomness 
(and hence attraction) leads to the formation of a bound state. The behavior of the 
bound state energy can then be used to extract an exponent of ( = 2/3 for superdiffusive 
wandering. By contrast, the replicated paths encountered in the directed wave problem 
(such as the two paths considered for eq. ( |145| )), although interacting, cannot form a 
bound state |HJ, as it is inconsistent with unitarity. This result also emerges in a natural 
way from our model of directed waves. In d = 2, for instance, it is easy to check that 
W2(r) ~ (1 — e5 rj o) _1 is the eigenstate of largest eigenvalue for the evolution of the relative 
coordinate. Hence, as t — > oo, for randomness 5-correlated in space and time, there is no 
bound state. This result holds in d > 2 and is not modified by short-range correlations in 



the randomness. The probability-conserving nature of eq. ( 147 ) is crucial in this regard 



105|| as it precludes a u5 D (r) attraction in eq. ( |150| ). Small perturbations that violate the 



conservation of probability lead to the formation of a bound state. In the language of the 
renormalization group, the scaling of directed waves in random media is governed by a 
fixed point that is unstable with respect to changes that do not preserve a strictly unitary 
evolution. 

Subsequently, a number of authors have obtained additional results from the random 



S'-matrix model. Following our work, Friedberg and Yu ||106|| calculated the leading terms 
in the scaling laws for the beam center in d > 2, and also the next-order corrections. The 
analytical results are in agreement with those presented above. Cule and Shapir [|107| 
extended the methods of this section to compute the higher moments of the probability 
distribution for directed waves in random media. If this probability distribution is multi- 



fractal, as claimed in ref. [ |1UU| , the higher moments should obey new scaling laws whose 



exponents are not simply related to those of the lower moments. Within the framework of 
the S'-matrix model, Cule and Shapir did not find evidence for multifractal scaling, while 
suggesting that certain aspects of the scaling behavior may be sensitive to details of the 
unitary time evolution. 
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P. Summing all Paths in High Dimensions 

In the next few sections we shall return to the non-random Ising model. The high 
temperature series can be approximately summed so as to reproduce mean-field (gaussian) 
behavior. This correspondence provides a better understanding of why such behavior is 
applicable in high dimensions, and also prepares the way for the exact summation of the 
series in two dimensions (next section). We shall then use these methods to look at two 
dimensional random Ising models. The high temperature series for the partition function 
of the non-random Ising model on a <i-dimensional hypercubic lattice is 

Z=J2e K ^(^ a ^ =2 N cosh dN KxS , (154) 

in} 

where S is the sum over all allowed graphs on the lattice, each weighted by r raised to 
the number of bonds in the graph, with r = t&nhK. The allowed graphs have an even 
number of bonds per site. The simplest graphs have the topology of a single closed loop. 
There are also graphs composed of disconnected closed loops. Keeping in mind cumulant 
expansions, we set 

S = sum over contribution of all graphs with one loop , (155) 

and introduce another sum, 

S' = exp (S) =1 + H + - (S) 2 + - (S) 3 + • • • 

1 ' 2 l ' 6 l ' (156) 

=1 + (1 loop graphs) + (2 loop graphs) + (3 loop graphs) + • • • 

Despite their similarities, the sums S and S' are not identical in that S' includes 
additional graphs where a particular bond contributes more than once. (In the original 
sum S, each lattice bond contributes a factor of 1 or r. After raising S to a power n, 
a particular bond may contribute up to r n .) In a similar approximation, we shall allow 
additional closed paths in S in which a particular bond is traversed more than once in 
completing the loop. Qualitatively, S is the partition function of a gas of self-avoiding 
polymer loops with a monomer fugacity of r. The self-avoiding constraint is left out in 
the partition function 5", which thus corresponds to a gas of phantom polymer loops. The 
corresponding free energy is 

In S' =sum over all closed random walks on the lattice x r length of walk 



=Ny^ — (number of closed walks of £ steps starting and ending at 0), 



(157) 



l 

52 



where extensivity is guaranteed since (up to boundary effects) the same loop can be started 
from any point on the lattice. The additional factor of l/£ accounts for the £ possible 
starting points for a loop of length £. To count the number of paths we introduce a set of 
N x N matrices, 

(i|W(£)|j) = number of walks from j to i in £ steps . (158) 

Then 

In S' - e 



£y<0|W(*)|0> • (159) 



N 

e 

Similarly, the spin-spin correlation function 

(a(0)a(r)) = i^a(0)a(r)n(l + ^^) , (160) 

{°~i} (ij) 

is related to the sum over all paths connecting the points and r on the lattice. In 
addition to the simple paths that directly connect the two points, there are disconnected 
graphs that contain additional closed loops. In the same approximation of allowing all 
intersections between paths, the partition function S' can be factored out of the numerator 
and denominator of eq. ( |160| ) , and 

<a(0Mr>«XV<r|W(*)|0> . (161) 

i 

As the walks are Markovian, their number can be calculated recursively. First note 
that any walk from to r in £ steps can be regarded as a walk from to some other point 
r' in £ — 1 steps, followed by a single step from r' to r. Summing over all possible locations 
of the intermediate point leads to 

(r\W{£)\Q) =^(r\W{l)\r') x (r'\W{£ - 1)|0) 

(162) 
= <r|2W(*-l)|0) , 

where the sum corresponds to the product of two matrices and we have defined T = W(l). 
The recursion process can be continued and 

W(£)=TW(£-l)=T 2 W(£-2) 2 = --- = T e . (163) 

Thus all lattice random walks are generated by the transfer matrix T, whose elements are 

/ it-" I '\ _ / 1 if r an d r' are nearest neighbors (1641 

1 otherwise 
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For example in d = 2, 



(x,y\T\x',y') = 5 y , y >(5 XjX ' +1 +8 X)X >-i) + 5 x , x >(5 y , y >+i + 6 y>y >- 1 ) 



(165) 



and successive actions of T on a walker starting at the origin, \x,y >= 6 Xt o5 yj o, generate 
the patterns 
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The value at each site is the number of walks ending at that point after £ steps. 

Various properties of random walks can be deduced from diagonalizing the matrix 
T. Due to the translational symmetry of the lattice, this is achieved in the Fourier basis 
(r|q) = e iqr /viV. For example, in d = 2 starting from eq (|165| ), it can be checked that 



(x 



y\T\q x , q y ) = Y^ ( x > y\ T \ x ', y') ( x ', y'lix, q v ) 



e „„ UW1 + ^M + e- e -*to + i. + ^(.-i. 



(166) 






i(qxX+q y y) 



[2cosg x + 2cos<? y ] =T(q xl q y ) (x,y\q x ,q y ) 



The generalized eigenvalue for a <i-dimensional hypercubic lattice is 

d 



T(q) = 2^2 cos q a . 

a=l 

The correlation function in eq. ( |161|) is now evaluated as 

oo oo 

(cr(r)cr(O)) « Y, T " < r I^WI°> = £ < r l^ T )'l°> 

£ £ 



(167) 




(168) 



iqr 



--N 



(27r) d N l-2rj:Li^sq a J (2yr) d 1 - 2r£ Q cosg, 



d d q 



3 iqr 



For r — > 0, the shortest path costs least energy and (cr(O)cr(r)) ~ r' r '. As r increases, 
larger paths dominate the sum because they are more numerous (i.e. entropically favored). 
Eventually there is a singularity for 1 — rT(0) = 0, i.e. at 2d x r c = 1, when arbitrarily 



54 



long paths become important. For r < r c , the partition function is dominated by small 
loops, and a polymer connecting two far away points is stretched by its line tension. When 
the fugacity exceeds r c , the line tension vanishes and loops of arbitrary size are generated. 
Clearly the neglect of intersections (which leads to a finite density) is no longer justified 
in this limit. This transition is the manifestation of Ising ordering in the language of 
paths representing the high temperature series. On approaching the transition from the 
high temperature side, the sums are dominated by very long paths. Accordingly, the 
denominator of eq. ( |168| ) can be expanded for small q as 

d 
1 - rT(q) = 1 - 2r ^ cos q a « (1 - 2dr) + rq 2 + 0(q 4 ) « t c (C 2 + q 2 + 0(q 4 )), (169) 

a=l 



where 

£=(— — ) .. (HO) 



l-2c^- 1/2 



The resulting correlation functions, (cr(O)cr(r)) oc f , 2 ? d e tq ' r /(g 2 + £~ 2 ), are identical to 
those obtained from a free (gaussian) field theory, and 

for r <C £ 



(a(0Mr))*< £* ■ (171) 



r (d-l)/2 



for r ^> £ 



The correlation length in eq. (|170|) diverges as ^ ~ (r c — t) ", with the exponent ^ = 1/2. 
We can also calculate the free energy in eq. ( |159| ) as 



00 J. 



UlS ' =J2 T T (o\w(i)\o) = (o 



N 



E T ' r 



e 



o 



= - (0| ln(l - rT)|0) = -N J -±± (0|q) In (1 - rT(q)) (q|0) (172) 

=-} (2^ ln ^-^E C0S ^J • 

In the vicinity of the critical point at r c = l/(2<i) the argument of the logarithm is propor- 
tional to (q 2 + £,~ 2 ) from eq. ( |169| ). This is precisely the free energy of a free field theory, 
and scales as 

/sing « r d cx (r c - r) d / 2 . (173) 
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The singular part of the heat capacity, obtained after taking two derivatives, is governed 
by the exponent a = 2 — d/2. Note that in evaluating the sums appearing in eqs. (|168| ) 



and ( |172|) the lower limit for £ is not treated carefully. The series in eq. (|168| ) is assumed 



to start from £ = 0, and that of eq. ( |172| ) from £ = 1. In fact the first few terms of both 



series may be zero because the number of steps is not sufficient to reach from r, or to 
from a closed loop. This is not a serious omission in that the singular behavior of a series 
is not effected by its first few terms. Treating the first few terms properly can only add 
analytic corrections to the singular forms calculated in eqs. (|168|) and ( |172| ). 

The equivalence of these results to a free field theory is a manifestation of field-particle 
duality. In a field theoretical description, (imaginary) time appears as an additional dimen- 
sion, and the two point correlations describe the probability of propagating a particle from 
one point in space-time to another. In a wave description, this probability is calculated by 
evolving the wave function using the Schrodinger equation. Alternatively, the probability 
can be calculated as the sum over all (Feynman) paths connecting the two points, each 
path weighted with the correct action. The second sum is similar to the above calculation 
of (a(r)a(O)). 

This approach provides an interesting geometrical interpretation of the phase tran- 
sition. The establishment of long range order implies that all parts of the system have 
selected the same state. This information is carried by the bonds connecting nearest 
neighbors, and can be passed from the origin to a point r through all paths connecting 
these two points. The fugacity r is a measure of the reliability of information transfer 
between neighboring sites. Along a one dimensional chain, unless r = 1, the transferred 
information decays at large distances and it is impossible to establish long range order. 
In higher dimensions there are many more paths, and by accumulating the information 
from all paths, it is possible to establish order at r c < 1. Since the number of paths of 
length £ grows as (2d) 1 while their information content decays as r , the transition occurs 
at r c = l/(2d). (A better approximation is obtained by including some of the constraints 
by noting that the random walk cannot back track. In this case the number of walks 
grows as (2d — 1)^.) The total information from paths of length £ is weighted by (2dr)^, 
and decays exponentially for r < r c . The characteristic path length, £ = — l/ln(2<ir), 
diverges as (r c — t) _1 on approaching the transition. For paths of size £ <C £ there is very 
good information transfer. Such paths execute random walks on the lattice and cover a 

- 1/2 

distance £ ~ £ . The divergence of 
the random walk nature of the paths. 



- 1/2 

distance £ ~ £ . The divergence of v with an exponent of 1/2 is thus a consequence of 
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Why does the classical picture fail for d < 4? Let us focus on the dominant paths close 
to the phase transition. Is it justified to ignore the intersections of such paths? Random 
walks can be regarded as geometrical entities of fractal (Hausdorf) dimension df = 2. 
This follows from the general definition of dimension relating the mass and extent of an 
object by M oc R df . The size of a random walk (R oc £) is the square root of its length 
(M oc £). Two geometrical entities of dimensions d± and d<2 will generally intersect in 
(i-dimensional space if d\ + d^ > d. Thus our random walkers are unlikely to intersect in 
d> d u = 2 + 2 = 4, and the results obtained by neglecting intersections are asymptotically 
valid. Below the upper critical dimension of 4, random walks have frequent encounters and 
their intersections must be treated correctly. The diagrams obtained in the perturbative 
calculation of the propagator in a ^> 4 theory correspond precisely to taking into account 
the intersections of paths. (Each vertex corresponds to one intersection.) It is now clear 
that the constraint of self-avoidance will swell the paths beyond their random walk size 
leading to an increase in the exponent v. Below the transition, the length of dominant 
paths grows without bound and the self-avoiding constraint is always necessary. 

Q. The Ising Model on a Square Lattice 



As indicated in eq. ( |154| ), the Ising partition function is related to a sum S, over 
collections of paths on the lattice. The allowed graphs for a square lattice have 2 or 4 
bonds per site. Each bond can appear only once in each graph, contributing a factor of 
r = tannic. While it is tempting to replace S with the exactly calculable sum S", of all 
loops of random walks on the lattice, this leads to an overestimation of S. The differences 
between the two sums arise from intersections of random walks and can be divided into 
two categories: 

(a) There is an over-counting of graphs which intersect at a site, i.e. with 4 bonds through 
a point. Consider a graph composed of two loops meeting at a site. Since a walker 
entering the intersection has three choices, this graph can be represented by three 
distinct random walks. One choice leads to two disconnected loops; the other two are 
single loops with or without a self-crossing in the walker's path. 

(b) The independent random walkers in S' may go through a particular lattice bond more 
than once. 

Including these constraints amounts to introducing interactions between paths. The 
resulting interacting random walkers are non-Markovian, as each step is no longer inde- 
pendent of previous ones and of other walkers. While such interacting walks are not in 
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general amenable to exact treatment, in two dimensions an interesting topological property 
allows us to make the following assertion: 

S = > collections of loops of random walks with no U turns 

^ (174) 

number of bonds / -in number of crossings 

The negative signs for some terms reduce the overestimate and render the exact sum. 
Proof: We shall deal in turn with the two problems mentioned above 
(a) Consider a graph with many intersections and focus on a particular one. A walker 
must enter and leave such an intersection twice. This can be done in three ways, only one 
of which involves the path of the walker crossing itself (when the walker proceeds straight 
through the intersection). This configuration carries an additional factor of (-1) according 
to eq. (|174j ). Thus, independent of other crossings, these three configurations sum up to 



contribute a factor of 1. By repeating this reasoning at each intersection we see that the 
over-counting problem is removed, and the sum over all possible ways of tracing the graph 
leads to the correct factor of one. 

(b) Consider a bond that is crossed by two walkers (or twice by the same walker) . We can 
imagine the bond as an avenue with two sides. For each configuration in which the two 
paths enter and leave on the same side of the avenue, there is another one in which the 
paths go to the opposite side. The latter involves a crossing of paths and hence carries a 
minus sign with respect to the former. The two possibilities thus cancel out! The reasoning 
can be generalized to multiple passes through any bond. The only exception is when the 
doubled bond is created as a result of a U-turn. This is why such backward steps are 
explicitly excluded from eq. (|174j) . 

Let us label random walkers with no U-turns, and weighted by (_i) number of crossings^ 
as RW*s. Then as in eq. ( p.56[ ) the terms in S can be organized as 

S = ^(RW*s with 1 loop) + ^(RW*s with 2 loops) + ^(RW*s with 3 loops) + • • • 

= exp i^(RW*s with 1 loop)! . 

(175) 
The exponentiation of the sum is justified since the only interaction between RW*s is the 
sign related to their crossings. As two RW* loops always cross an even number of times, 
this is equivalent to no interaction at all. Using eq. (|154|) , the full Ising free energy is 
calculated as 

In Z = N In 2 + 2N In cosh K + J^ (RW*s with 1 loop x t* of bonds ) . (176) 
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Organizing the sum in terms of the number of bonds, and taking advantage of the trans- 
lational symmetry of the lattice (up to corrections due to boundaries), 

-^ = ln(2cosh 2 K)+^^(0|W*(£)|0) , (177) 

where 

(0|W*(^)|0) =number of closed loops of I steps, with no U turns, from to 

(!78) 

/ 1\# °f crossings 

The absence of U-turns, a local constraint, does not complicate the counting of walks. 
On the other hand, the number of crossings depends on all sites previously crossed by the 
walker and is a non-Markovian property. Fortunately, in two dimensions it is possible to 
obtain the parity of the number of crossings from local considerations. The first step is to 
make the walks oriented by placing an arrow along the direction that the path is traversed. 
Since any path can be traversed in two directions, 

(0\W*(£)\0) = - ^2 oriented RW* loops of £ steps, no U turns, from to x (-l) nc , 

(179) 
where n c is the number of self-crossings of the loop. We can now take advantage of the 
following topological result ||108|| : 

Whitney's Theorem: The number of self-crossings of a planar loop is related to the total 
angle O, through which the tangent vector turns in going around the loop by 

K) mo d 2 = ( 1 + ^- J • (180) 

V 27 Vmod2 

This theorem can be checked by a few examples. A single loop corresponds to O = ±27r, 
while a single intersection results in O = 0. 

Since the total angle O is the sum of the angles through which the walker turns at 
each step, the parity of crossings can be obtained using local information alone, as 

(_l) n « = e *™c = exp itt (l + — J = -e^S=^ J , (181) 

where 6j is the angle through which the walker turns on the j th step. Hence 

(0|W*(^)|0) = — 2_, oriented RW* loops of £ steps, with no U turns, from to 

x exp I - 2_. local change of angle by the tangent vector J 

(182) 
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The angle turned can be calculated at each site if we keep track of the directions of arrival 
and departure of the path. To this end, we introduce a label \x for the 4 directions going 
out of each site, e.g. \i = 1 for right, \x = 2 for up, \i = 3 for left, and \x = 4 for down. We 
next introduce a set of 4iV x 4iV matrices generalizing eq. ( |158| ) to 

(x 2 y 2 , ftelW* (£)\xiyi, Hi) = / oriented random walks of £ steps, with no U turns, 

departing (xi,yi) along hi, proceeding along [i 2 after reaching (#2,2/2) x e2 i=1 J 

(183) 
Thus /X2 specifies a direction taken after the walker reaches its destination. It serves to 
exclude some paths (e.g., arriving along -fJ, 2 ), and leads to an additional phase. As in 
eq. (|162|) , due to their Markovian property, these matrices can be calculated recursively as 

(x 2 y2i H2\W* (tylxxyx, (ii) = 
^2 (x 2 y 2l fj> 2 \T*\x'y\ij') (x'y',fj,'\W*(£- l)\xiyi, m) = ( 184 ) 

x'y' ,fj,' 

(x 2 y 2 ,n 2 \T*W*(£- l)|a;ij/i,/*i> = (x 2 y2,fJ>2\T* e \x 1 yi, //1) , 

where T* = W*(l) describes one step of the walk. The direction of arrival uniquely 
determines the nearest neighbor from which the walker departed, and the angle between 
the two directions fixes the phase of the matrix element. We thus generalize eq. ( [165) to 



a 4 x 4 matrix that keeps track of both connectivity and phase between pairs of sites, i.e. 

(x'y'\T*\xy) = 

(x'y'\x + ly) (x'y'\x + ly) e^ (x'y'\x + ly) e _i ? 

(x'y'\xy + 1) e - ^ (x'y'\xy + l) (x'y'\xy + 1) e^ 

(x'y'\x — ly) e~~ (x'y'\x — ly) (x'y'\x — ly) e~ 

(x'y'\xy — 1) ei (x'y'\xy — 1) e _i ? (x'y'\xy — 1) 

(185) 

where < xy\x'y' >= b~ x>x '0~ y ,y' ■ 

Because of its translational symmetry, the AN x AN matrix takes a block diagonal 
form in the Fourier basis (xy\q x q y ) = e l ^ qxX+qyy ' / 'yN , i.e. 

^2(x'y',H'\T*\xy,fi)(xy\q x q y ) = (//|T*(q)|/i) (x'y'\q x q y ) . (186) 

xy 

Each 4x4 block is labelled by a wavevector q = (q x , q y ), and takes the form 

T *( C l) = n J(q x ~^) „iq x „i(q* + ?r) " ( 187 ) 
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e -i(«a,-f ) 





e -*(9« + f)-i 


;-*(«» + f) 


e -*9y 


e -*(«l/-f ) 








g*(fc-f ) 


gi«2x 


e %*+i) 


e *fov + f) 





e *(9«-f) 


e *5y 



To ensure that a path that starts at the origin completes a loop properly, the final 
arrival direction at the origin must coincide with the original one. Summing over all 4 such 
directions, the total number of such loops is obtained from 

4 1 1 

(0\W*(£)\0) = J2 (00,/*|T**|00,m) = T7 E {xy^\T* l \xy,n) = ^tr(T^) . (188) 



ft=i 



xy,n 



Using eq. (|177| ), the free energy is calculated as 

'" Z : ln(2 cosh 2 K) - ^E^-(0|^*(£)|0) = ln(2cosh 2 i^)--^-tr 



N 



E 



rp*e T e 



1 



ln(2 cosh 2 K) + — tr ln(l - tT*) 
ln(2cosh 2 K) + -l-Etrln(l-rT*(q)) 



(189) 



But for any matrix M, with eigenvalues {A a }, 



trlnM = ^lnA a = lnJ]A a = IndetM 



Converting the sum over q in eq. ( |189| ) to an integral leads to 
InZ 



N 



ln(2cosh 2 J FQ + 



1 f d 2 q 



2 J (2tt) 2 



In < det 



1 _ re~ iqx —re-^g^-f ) o _ re -*(fe+i ) 

_ re -*(3s/ + f) l _ Te - i( ly _ re -*(9«-f) 




. re *(<iv+f) 



-re 







re 



iq.T 



-re 



i(9x + f) 



-re 



*(«f»-f) i 



re 



"/;/ 



Evaluation of the above determinant is straightforward, and the final result is 
InZ . ,_ . , „, If d 2 q 



N 



ln(2 cosh 2 K) + ~ / -^ In [(1 + r 2 ) 2 - 2r (1 - r 2 ) (cos q x + cos ^ 



(190) 



(191) 



Taking advantage of trigonometric identities, the result can be simplified to 
lHZ ln2+ * f^ ^^ln[co S h 2 (2K)- S mh(2K)(co S q x + co S q y )} 



N 



U-TT (27T) 



(192) 



While it is possible to obtain a closed form expression by performing the integrals exactly, 
the final expression involves a hypergeometric function and is not any more illuminating. 
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R. Singular Behavior 

To uncover the singularity in the free energy of the two dimensional Ising model in 
eq. (|191| ), we start with the simpler expression obtained by the unrestricted sum over 
random walks in eq. ( |172| ), (specializing to d — 2) 

/ G = ln(2cosh 2 K)- /"^^ln[l-2r(cos^ + cos^)] . (193) 

Apart from the argument of the logarithm, this expression is similar to the exact result. 
Is it possible that such similar functional forms lead to distinct singular behaviors? The 
singularity results from the vanishing of the argument of the logarithm at r c = 1/4. In the 
vicinity of this point we make an expansion as in eq. (|169| ), 

A g (t, q) = (1 - 4r) + rq 2 + 0(q 4 ) « r c (q 2 + Ay\ , (194) 

where Sr = t c — r. The singular part of eq. ( |193| ) can be obtained by focusing on the 
behavior of the integrand as q — > 0, and replacing the square Brillouin zone for the range 
of the integral with a circle of radius A rj 2-7T, 



i 



/sing. =- / -r^ln(g 2 +4-) 



, 2 + 4 ^ln^ 2 + 4 ' r/Tc 



4n \* r c J \ e 

Only the expression evaluated at q = is singular, and 

1 f St\ , {St 



(195) 



J o 



/-* = "UJ ln U' (196) 



The resulting heat capacity Cq oc d 2 fo/d 2 t, diverges as 1/St. Since eq. ( p.93| ) is not valid 
for r > r c , we cannot obtain the behavior of heat capacity on the low temperature side. 
For the exact result of eq. ( |191| ) , the argument of the logarithm is 

A*(r,q) = (l + r 2 ) 2 -2r(l-r 2 )(cosQ x + cos^) . (197) 

The minimum value of this expression, for q = 0, is 

A*(t,0) = (1 + t 2 ) 2 -4t(1-t 2 ) = (1-t 2 ) 2 +4t 2 -4t(1-t 2 ) = (l-r 2 -2r) 2 . (198) 
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Since this expression (and hence the argument of the logarithm) is always non-negative, 

the exact result is valid 



the integral exists for all values of r. As required, unlike eq. ([19c 

at all temperatures. There is a singularity when the argument vanishes for 

t 2 + 2t c -1 = =► r c = -l± v / 2 . (199) 

The positive solution describes a ferromagnet and leads to a value of K c = ln(v2 + l)/2. 
Setting St = r — r c , and expanding eq. (|197Q in the vicinity of q — > gives 



A*(r,q)^[(-2r c -2)5r] 2 + r c (l 



rlW + 



'iTi 



^O 



(200) 



The important difference from eq. ( [I94[ ) is that (St/t c ) appears at quadratic order. Fol- 
lowing the steps in eqs. (|195| ) and ( |196| ) , the singular part of the free energy is 



InZ 



N 



sing. 



A 



2ixqdq 
4n 2 



In 



q 2 +4 



9' 
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( 9 {St 




9 + 4 ~ 


87T 


V \Tc 


i( 


H)\ 


5r 


TV V 


rc) 


T c 



In 



q 2 + 4(5t/t c 



(201) 



The heat capacity is obtained by taking two derivatives and diverges as C(5r) s i ng . = 
A±ln|dr|. The logarithmic singularity corresponds to the limit a = 0; the peak is sym- 
metric and characterized by the amplitude ratio of A + /A- = 1. 

The graphical method presented in this section was originally developed by Kac and 
Ward ||109|| . The main ingredient of the derivation is the result that the correct accounting 
of the paths can be achieved by including a factor of (—1) for each intersection. (This was 
originally a conjecture by Feynman ||110|| , later proved by Sherman ||111|| .) The change of 
sign is reminiscent of the exchange factor between fermions, and indeed the final result can 
be obtained by mappings to free fermions ||112|| . 

In addition to the partition function, the correlation functions < o-jOj > can also be 
calculated by summing over paths [|113| . Since the combination q 2 +4(5r/r c ) 2 in eq. (|200| ) 
describes the behavior of these random walks, we expect a correlation length £ ~ \t c /St\, 
i.e. diverging with an exponent v — 1 on both sides of the phase transition, with an 
amplitude ratio of unity. The exponents a and v are related by the hyperscaling identity 
a = 2 — 2v. The critical correlations at r c are more subtle and decay as < o\Oj > c ~ 
l/|i — jl 77 with r\ = 1/4. Integrating the correlation functions yields the susceptibility, 
which diverges as x± — C±|^ r l~ 7 ; with 7 = 7/4 and C + /C- = 1. 
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S. The Two Dimensional Spin Glass 

The key to the exact solution of the non-random two dimensional Ising model pre- 
sented in the previous section is the reduction of the graphical expansion to sums over 
non-interacting (fermionic) random walks. This reduction depends on the geometrical 
properties of paths in d = 2 and is independent of the uniformity of the bonds J^. 
The same method can be applied to the random bond problem, reducing the problem 
to (fermionic) random walks in a random medium. This is an undirected version of the 
problem of DPRM, extensively mentioned in these notes. Due to randomness, we can no 
longer diagonalize the transfer matrix that generates these walks by Fourier transforma- 
tion. However, we can still examine such walks numerically by successive multiplications 
of the transfer matrix. Thus all random bond Ising models in d = 2 can in principle be 
solved exactly in polynomial time in their size L. This is not true for the three dimensional 
versions. 

We shall demonstrate the potentials of such exact methods by developing an algo- 
rithm for calculating the partition function of the ±J spin glass [|51 |. This is the model 



introduced in sec. I, where its high temperature correlation functions were examined with 
directed paths. The advantage of the model (or any other random mixture of + J, —J, and 
absent bonds) is that all computations can be performed in integer form, thus avoiding any 
floating point errors. Of course the main interest in the spin glass problem stems from the 
complexity of its low temperature states. Despite the great deal of work on spin glasses 



over the past decades [4,114 — 116] , the description of the phase transition and the nature 



of the glassy state remain controversial subjects |117| , |118| . Interactions with infinite-range 



|118(] lead to a solution with broken replica symmetry |4[ . It is not known, however, to what 
extent this mean-field result captures the behavior of short-range interactions ||118| , |120f . 
Monte Carlo simulations of spin glass problems are notoriously difficult due the ease with 
which the system can get trapped in long-lived metastable states. It is thus quite valuable 
to provide some exact information about the equilibrium low temperature behavior of spin 
glasses, even if that knowledge is limited to two dimensions. 



We start with the Edwards- Anderson (EA) Hamiltonian ||121 



H = ^J lJ a l a 3 , (202) 

(ij) 
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where the nearest neighbor quenched random bonds Jy are chosen from the bimodal 
distribution 

p{Jij) = \KJij-J) + \KJii + J) , (203) 

with J > 0. On a lattice with periodic boundary conditions (BCs), there are exactly 2N 
bonds, with N = L 2 the total number of spins. The high temperature expansion for the 
partition function takes the form 

27V 

Z = 2 N cosh 2N K Y^AiT 1 , (204) 

1=0 

where the coefficients A\ are pure integers. Note that Ag = for odd values of £ since 
closed loops on the square lattice necessarily traverse an even number of bonds. 

We can use the diagrammatic method introduced in the previous sections to transform 
the problem of summing the high temperature series into one of evaluating a local random 
walk. Every step proceeds exactly as before up to eq. (|162j ). However, the 4N x 4iV 



transfer matrix in eq. (|185| ) has to be modified to take care of the randomness in bonds. 
Each element of the matrix that connects sites i and j has to be multiplied with the 
reduced bond variable Sij = Sji = Jij/J, equal to +1 for ferromagnetic bonds and —1 for 
antiferromagnetic bonds. There is an additional complication in dealing with finite lattice 
size. For simplicity we shall use periodic boundary conditions on an L x L lattice. We 
must then take proper account of diagrams which wrap around the lattice. The correct 



result, based on the combinatorics of closed loops on periodic lattices ||122|| is Z = {—Z\ + 
Z 2 + Z 3 + Z 4 )/2 with 

Z x = 2 N cosh 2Ar K^/det [1 - T*t] . (205) 

Here, T* is the original 4N x AN transfer matrix, while T% , Tg, (and T£) are obtained 
respectively by changing the sign of a horizontal, vertical, (or both) column of bonds. 
The linear combination {—Z\ + Z 2 + Z% + Z^)j2 ensures that all diagrams, including those 
which loop the entire lattice, are weighted correctly in the final expression for the partition 
function. 

We have implemented ||123fl this algorithm on the computer as follows. Given a set of 



bonds {Jij}, we first construct the 4N x AN matrices T£ and compute the traces ti(T£ ) for 
£ < N. This step of the algorithm is the most computationally intensive. The coefficients 
of the series expansions for In Z\ are related to the traces by 



InZ 



1 °° 1 
ln[2cosh 2 K]-— £ 7 tr(T*V • (206) 



N l J 2N 
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Next, we compute the high temperature series for Z. This is done by exponentiating 
the series for lnZ\, followed by taking the linear combination that incorporates periodic 
boundary conditions. The high temperature expansion for Z is a polynomial in r with 
integer coefficients; the last term, of order 2N, is derived from the graph that traverses 
every bond on the square lattice. These 2N coefficients have an end-to-end symmetry 
that enables one to compute them from the first N powers of the transfer matrix. Finally, 
we expand powers of coshET and tannic and rewrite Z as a polynomial in e~^ J ; the end 
result Z = '^2 E g(E)e~ f3E yields the density of states. For an Ising model with ±J bonds, 
we can perform all these operations using only integer arithmetic. 

The algorithm has several desirable features. First, it returns the partition function Z 
as an exact integer result. In this way, it not only avoids the statistical uncertainties 
inherent in Monte Carlo simulation; it also avoids the floating point errors that creep into 
numerically "exact" calculations of Z in large systems at low temperatures. The algorithm 
thus provides us with an efficient and reliable way to investigate energies and entropies 
at low temperatures. This is particularly important in a system that exhibits a phase 
transition at T = 0, such as the ±J spin glass. We can also calculate other quantities, 
such as the roots of the partition function in the complex plane, or the number of low- 
level excitations, that are otherwise inaccessible. Unfortunately, the necessity of handling 
large integers (of order 2^) complicates what would otherwise be a rather straightforward 
algorithm to implement. 

A second advantage of the algorithm is that it executes in polynomial time. We 
estimate the algorithm's performance as follows. Computing the traces requires 0(N 3 ) 
arithmetic operations on integers of order 2 , while in general, the power series manipula- 
tions take much less time. We therefore expect the computation time to scale as r ~ N s , 
with 3 < S < 4. This stands in contrast to the numerical column to column transfer 



matrix method of Morgenstern and Binder ||124j| which has time and memory requirements 
that grow exponentially with system size. We obtained most of our results on dedicated 
DEC 3100 workstations. The largest lattice that we examined had N = 36 x 36 spins. 
Finally, we mention that the computation of the traces can be broken down into O(N) 
independent computations, so that a faster, parallel implementation of the algorithm (on 
a supercomputer or spread across several workstations) should be possible. 
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T. Results for the Two Dimensional Spin Glass 

We examined the ± J spin glass on lattices of size L = 4 to L = 36. Several realizations 
of randomness were studied: 8000 for L = 4, 6, 8; 2000 for L = 10, 12, 14; 800 for L = 16, 
18; 80 for L = 20, 22, 24; and 4 for L = 32, 36. We performed quenched averages by 
assigning an equal probability to each random sample: 9 = (1/S) J2 s @ s - To re duce the 
amount of statistical error, we only considered lattices in which exactly half the plaquettes 



were frustrated [124]. We also found it quite revealing to compare our results with those 



for the fully frustrated Ising model [|125| , |126|1 , as both models undergo phase transitions at 



T = 0. The typical output of our algorithm is a set of integers g(E) for the number of 
states. Using the density of states we can perform calculations in either the microcanonical 
or canonical ensemble. In the limit of infinite size, of course, the two ensembles should 
yield identical results. 

We used the algorithm first to study the thermodynamic properties of the ground 
state, fitting the data for the ground-state energy and entropy to the finite-size form 
II = foo + a/L 2 . Extrapolating to infinite lattice size, we estimate eo/J = —1.404 ± 0.002 
and so = 0.075 ± 0.002. These results are consistent with previous MC [|127| , |128|| and 



column to column transfer matrix ||129|| estimates. We also used the algorithm to study 
the number of low-level excitations. On a lattice with periodic BCs, the lowest excited 
state has an energy AJ above the ground state. The quantity e AS = g(Eo + 4J)/g(Eo) 



measures the degeneracy ratio of these excited states. We find that ASsg grows faster 
than IniV indicating that the low-lying excitations of the ±J spin glass involve spin flips 
on large length scales. 

The abundance of low-lying excitations affects the low-temperature behavior of the 
heat capacity. In a finite system with energy gap 4 J, the heat capacity vanishes as C ~ 
(3 2 e~ 4l3J . As pointed out by Wang and Swendsen [ 127|| , this behavior can break down 



in the thermodynamic limit. The ID Ising model with periodic BCs shows how this can 
happen: the energy gap is 4J, but the heat capacity of an infinite system vanishes as 
C\d ~ (3 2 e~ 2 @ J . The anomalous exponent reflects the fact that the number of lowest 
excited states grows as iV 2 . From Monte Carlo and column to column transfer matrix 
studies, Wang and Swendsen |127|] conclude that Csg ~ (3 2 e~ 2f3J for the 2D ±J spin 



glass as well. Our results find a disagreement in slope between ASm and ASsg versus 
N, leading us to suggest a different form for Csg- As motivation, we appeal to another 
exactly soluble model with a phase transition at T = 0: the fully frustrated (FF) Ising 
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model on a square lattice [125]. On a periodic lattice, the lowest excited states of the FF 
model have energy 4 J above the ground state. The large number of low-lying excitations, 
however, causes the heat capacity to vanish as Cff ~ /3 3 e _4 ^ J . Note the extra power of 
temperature. Our data compare much better to ASff than to ASid, suggesting a similar 
behavior may describe the ± J spin glass, e.g. Csg ~ (3 2+p e~ 4:,3J with p ^ 0. As we shall 
see below, there are other reasons to favor this form. 

One way to investigate phase transitions is to examine the roots of the partition 



function Z in the complex plane. This was first done by Fisher [|130|1 in a study of the 



2D Ising model with uniform + J bonds. Fisher calculated the distribution of roots of 
the partition function in the complex z = e~ 2,3J plane. In the limit of infinite lattice 
size, he showed that the roots condense onto two circles centered at z = ±1, and related 
the singular behavior in the free energy to the distribution of roots in the vicinity of the 
positive real axis. Since a system of finite size does not exhibit non-analytic behavior, it is 
clear that the roots of the partition function can only close in on the positive real axis in 
the thermodynamic limit. The zeros of partition functions are thus subject to finite-size 
scaling ||131|| . At a finite-temperature phase transition, the complex zero T(L) closest to 



the positive real axis obeys \T(L) — T c \ ~ L~ yt ; likewise, the correlation length diverges 
as £ ~ (T — T c )~ v with v = 1/yt- On the other hand, at a T = phase transition, such 
as occurs in the ID Ising model, one finds |z(£)| ~ L~ Vz with £ ~ z~ 1 ' Vz 1 where z(L) is 
the smallest root in the complex z = e~ 2 @ J plane. For example, the partition function of 
a ID Ising chain with periodic boundary conditions in eq. ( |28| ) , 

Z = 2 L cosh L (/3J) [l + tanh L (/3J)j . (207) 

has its smallest root at tanh(/3J) = e ±l7r ' L , or z(L) = ±z tan(7r/2L). As L — >• oo, the 
magnitude of the root scales as z(L) ~ L~ , consistent with the fact that the correlation 
length in the ID Ising model diverges as £ ~ e 2/3J . 

In the ±J spin glass, we observed that, for most realizations of randomness, the 
smallest root z(L) falls on the imaginary axis. One might expect that the probability 
distribution for the magnitude of this root assumes a scale-invariant form as L — > oo. We 
were unable to verify this hypothesis due to insufficient data on large lattices. Instead, 
we examined the statistics of w(L), where u = z 2 = e~ 4 @ J . On a square lattice with 
periodic boundary conditions, the partition function for a ±J spin glass is polynomial in 
e -4/3j ^Y e therefore looked at the scaling of roots in the complex u = e~ 4)3J plane. The 
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results could be fitted to |w(X)| ~ L~ 2 - 2 with y u = 2.2 ± 0.1; this suggests to us that 
the correlation length in the ±J spin glass diverges as £ ~ e 2 @ J . Additional powers of 
temperature and/or finite-size effects might explain the slight deviation from y u = 2. Note 
that this behavior for the correlation length is consistent with the hyperscaling relation 
lnZ S i n g U i ar oc £ -2 , and our claim that, up to powers of temperature, the heat capacity 
diverges as C ~ e~ 4,3J . This result disagrees with previous studies ||127| , |129j , |132| that 
report £ ~ T~ u ', with v « 2.6 — 2.8. An analytical approach, closely related to averaging 



the fermionic path integrals ||133|| , does support exponential correlations, but with a gap 
of 4 J. 

A great deal of information on spin glasses has been obtained by examining 'defects' 
(droplets) in finite systems. The cost of a defect of size L is related to the difference in 
free energies with periodic and anti-periodic BCs. At T = 0, this reduces to the difference 
in energy between the ground states. On an L x L lattice, the defect energy measures the 
effective block coupling p4] , |iT?yTTgl , |T2ll , p51 J' on length scale L. Let p(L) be the fraction 
of L x L blocks for which J' ^ 0. Scaling arguments [ |135| suggest that p(L) ~ L~ v , where 
r\ is the critical exponent that characterizes the power law decay of correlations < uqUl > 2 
at T = 0. Plotting p(L) versus L, we find r\ = 0.22 ± 0.06 in agreement with previous 



results [ p.35[ , |127|| . Besides the defect energy, we also looked at the defect entropy SSl, i-e. 



the difference in zero-temperature entropies with periodic and anti-periodic BCs. The data 



could be fitted to 5S 2 L ~ L 2ys with y s = 0.49 ± 0.02. This is curiously close to the result 
5Sl ~ L 1 / 2 , expected if entropy changes due to reversing the different bonds along the 
boundary are statistically independent. The defect entropy in the fully frustrated Ising 
model approaches a constant value with 1/L 2 corrections. It is straightforward, moreover, 
to show that in the ID Ising model the defect entropy scales as SSl ~ InL. Both these 
behaviors are markedly different from the spin glass. We do not know any obvious relation 
between the finite-size scaling of the defect entropy and other quantities at T = 0. More 
details of the algorithm and results can be found in ref. ||123| . 
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